ParallelDomainBase.cc 90.1 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
#include <algorithm>
2
3
#include <iostream>
#include <fstream>
Thomas Witkowski's avatar
Thomas Witkowski committed
4

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

#include "petscksp.h"
27
28
29

namespace AMDiS {

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

    return 0;
  }

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

43
  ParallelDomainBase::ParallelDomainBase(ProblemIterationInterface *iIF,
44
45
					 ProblemTimeInterface *tIF,
					 FiniteElemSpace *fe,
46
					 RefinementManager *refinementManager)
47
48
    : iterationIF(iIF),
      timeIF(tIF),
49
      name(iIF->getName()),
50
51
      feSpace(fe),
      mesh(fe->getMesh()),
52
      refineManager(refinementManager),
53
      initialPartitionMesh(true),
54
55
      d_nnz(NULL),
      o_nnz(NULL),
56
      nRankDofs(0),
57
      rstart(0),
58
      nComponents(1),
59
60
      deserialized(false),
      lastMeshChangeIndex(0)
61
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
62
63
64
65
66
67
68
    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");

69
70
71
72
73
74
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
  }

75
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
76
  {
77
78
79
80
81
82
83
    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)
84
85
      return;

86
87
88
89
90
    // 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();

91
92
    MSG("START PARTITIONING!\n");

93
94
95
96
97
98
99
    // 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
100
101
#if (DEBUG != 0)
    ElementIdxToDofs elMap;
102
    dbgCreateElementMap(elMap);
103
104
105
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
106

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

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

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

118
119
120
121
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
122
123
124
125
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

126
127
128
129
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

130
131
    // === If in debug mode, make some tests. ===

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

139
140
    debug::writeMesh(feSpace, -1, "macromesh");   
#endif
141

142
143
144
    // === Create periodic dof mapping, if there are periodic boundaries. ===

    createPeriodicMap();
145

146
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
147

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

Thomas Witkowski's avatar
Thomas Witkowski committed
151
    if (globalRefinement > 0) {
152
      refineManager->globalRefine(mesh, globalRefinement);
153

154
      updateLocalGlobalNumbering();
155

156
      // === Update periodic mapping, if there are periodic boundaries. ===
157

158
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
159
    }
160
161
  }

162
163

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
164
  {}
165

166
167
168
  
  void ParallelDomainBase::updateDofAdmins()
  {
169
170
171
    FUNCNAME("ParallelDomainBase::updateDofAdmins()");

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
172
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
173
174
175

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

184
185
186
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
187
188
189
    }
  }

190

191
192
193
194
195
196
197
198
199
200
  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)
201
	("Mesh is already refined! This does not work with parallelization!\n");
202
203
204
205
206
207
208
209
210
211
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

212

213
214
  void ParallelDomainBase::setDofMatrix(DOFMatrix* mat, int dispMult, 
					int dispAddRow, int dispAddCol)
215
  {
216
217
218
219
    FUNCNAME("ParallelDomainBase::setDofMatrix()");

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

220
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
221
222
223
224
225
226
    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());

227
    typedef traits::range_generator<row, Matrix>::type cursor_type;
228
229
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

230
231
232
233
234
    std::vector<int> cols;
    std::vector<double> values;
    cols.reserve(300);
    values.reserve(300);

235
236
237
    // === Traverse all rows of the dof matrix and insert row wise the values ===
    // === to the petsc matrix.                                               ===

238
239
240
241
242
243
    for (cursor_type cursor = begin<row>(mat->getBaseMatrix()), 
	   cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {

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

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

249
250
251

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

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

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

264
265
	  // 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.
266
 	  if (periodicRow == false && periodicDof.count(globalColDof) > 0) {
267
268
269
270
271
	    // 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.
272
 	    cols.push_back(colIndex);
273
 	    values.push_back(value(*icursor) * scalFactor);
274

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

290
291
292
293
294

      // === 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.
295
      int rowIndex = globalRowDof * dispMult + dispAddRow;
296
      
297
      if (periodicRow) {
298
299
300
	// The row dof is periodic, so send dof to all the corresponding rows.

	double scalFactor = 1.0 / (periodicDof[globalRowDof].size() + 1.0);
301
	
302
	int diagIndex = -1;
303
	for (int i = 0; i < static_cast<int>(values.size()); i++) {
304
305
306
307
308
	  // 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;
309
310
311
	  else
	    diagIndex = i;
	}
312
313
314
315
316
317
318
	
	// 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.
319
320
321
	if (diagIndex != -1)
	  values[diagIndex] = 0.0;

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

      } else {
331
332
333
	// 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);
334
      }    
335
    }
336
  }
337

338

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

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

	double value = *dofIt / (periodicDof[globalRow].size() + 1.0);
354
355
356
357
358
359
360
	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);
	}
361

362
      } else {
363
	// The dof index is not periodic.
364
365
366
	double value = *dofIt;
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
      }
367
    }    
368
369
  }

370

371
  void ParallelDomainBase::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
372
  {
373
    FUNCNAME("ParallelDomainBase::createPetscNnzStructure()");
374

375
376
    TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n");
    TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n");
377

378
379
380
381
382
383
    d_nnz = new int[nRankRows];
    o_nnz = new int[nRankRows];
    for (int i = 0; i < nRankRows; i++) {
      d_nnz[i] = 0;
      o_nnz[i] = 0;
    }
384

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

390
391
392
393
394
395
    // 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;

396
397
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
398
399
400
401
402
403
 	if ((*mat)[i][j]) {
	  Matrix bmat = (*mat)[i][j]->getBaseMatrix();

	  traits::col<Matrix>::type col(bmat);
	  traits::const_value<Matrix>::type value(bmat);
	  
404
	  typedef traits::range_generator<row, Matrix>::type cursor_type;
405
406
	  typedef traits::range_generator<nz, cursor_type>::type icursor_type;
	  
407
408
409
	  for (cursor_type cursor = begin<row>(bmat), 
		 cend = end<row>(bmat); cursor != cend; ++cursor) {

410
411
412
	    // 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;
413

414
	    if (isRankDof[*cursor]) {
415
416
417
418
419
420

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

#if (DEBUG != 0)    
423
	      if (localPetscRowIdx < 0 || localPetscRowIdx >= nRankRows) {
424
		std::cout << "ERROR in rank: " << mpiRank << std::endl;
425
426
		std::cout << "  Wrong r = " << localPetscRowIdx << " " << *cursor 
			  << " " << mapLocalGlobalDofs[*cursor] << " " 
427
428
429
430
			  << nRankRows << std::endl;
		ERROR_EXIT("Should not happen!\n");
	      }
#endif
431
	      
432
	      // Traverse all non zero entries in this row.
433
434
435
	      for (icursor_type icursor = begin<nz>(cursor), 
		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
		if (value(*icursor) != 0.0) {
436
		  int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
437

438
439
440
441
442
		  // 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]++;
443
		  else
444
		    o_nnz[localPetscRowIdx]++;
445
446
447
		}    
	      }
	    } else {
448
449
450
451
452
	      
	      // === 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.                          ===
453

454
455
	      // Find out who is the member of this DOF.
	      int sendToRank = -1;
456
457
458
459
460
	      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) {
461
462
463
464
465
466

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

467
468
469
470
471
472
473
474
475
476
477
		    sendToRank = it->first;
		    break;
		  }
		}

		if (sendToRank != -1)
		  break;
	      }

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

478
	      // Send all non zero entries to the member of the row DOF.
479
480
481
	      for (icursor_type icursor = begin<nz>(cursor), 
		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
		if (value(*icursor) != 0.0) {
482
		  int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
483
		  
484
485
		  sendMatrixEntry[sendToRank].
		    push_back(std::make_pair(petscRowIdx, petscColIdx));
486
487
488
		}
	      }

489
490
491
492
	    } // if (isRankDof[*cursor]) ... else ...
	  } // for each row in mat[i][j]
	} // if mat[i][j]
      } 
493
494
    }

495
    // === Send and recv the nnz row structure to/from other ranks. ===
496

Thomas Witkowski's avatar
Thomas Witkowski committed
497
    StdMpi<MatrixNnzEntry> stdMpi(mpiComm, true);
498
499
500
    stdMpi.send(sendMatrixEntry);
    stdMpi.recv(sendDofs);
    stdMpi.startCommunication<int>(MPI_INT);
Thomas Witkowski's avatar
Thomas Witkowski committed
501

502
503
    // === Evaluate the nnz structure this rank got from other ranks and add it to ===
    // === the PETSc nnz data structure.                                           ===
504

505
506
507
508
509
510
    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;
511

512
	  int localRowIdx = r - rstart * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
513

514
515
516
	  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);
517
	  
518
	  if (c < rstart * nComponents || c >= rstart * nComponents + nRankRows)
519
	    o_nnz[localRowIdx]++;
520
	  else
521
	    d_nnz[localRowIdx]++;
522
523
	}
      }
524
    }  
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
  }


  void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
  {
    FUNCNAME("ParallelDomainBase::fillPetscMatrix()");

    clock_t first = clock();

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

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

    if (!d_nnz)
      createPetscNnzStructure(mat);
Thomas Witkowski's avatar
Thomas Witkowski committed
550

551
    // === Create PETSc matrix with the computed nnz data structure. ===
552
553

    MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
554
555
		    0, d_nnz, 0, o_nnz, &petscMatrix);

556
557
    INFO(info, 8)("Fill petsc matrix 1 needed %.5f seconds\n", TIME_USED(first, clock()));

558
559
560
561
562
563
#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
564

565
    // === Transfer values from DOF matrices to the PETSc matrix. === 
566
567
568

    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
569
	if ((*mat)[i][j])
570
571
	  setDofMatrix((*mat)[i][j], nComponents, i, j);

572
573
    INFO(info, 8)("Fill petsc matrix 2 needed %.5f seconds\n", TIME_USED(first, clock()));

574
575
576
    MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);

577
578
    // === Transfer values from DOF vector to the PETSc vector. === 

579
    for (int i = 0; i < nComponents; i++)
580
      setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i);
581

582
583
584
    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);

585
    INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", TIME_USED(first, clock()));
586
587
588
  }


589
  void ParallelDomainBase::solvePetscMatrix(SystemVector &vec)
590
591
592
  {
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

593
594
595
596
597
598
599
600
601
#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

602
    // === Init Petsc solver. ===
603

604
    KSP solver;
605
606
    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 
607
    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
608
609
    KSPSetType(solver, KSPBCGS);
    KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0);
610
    KSPSetFromOptions(solver);
611
612
    // Do not delete the solution vector, use it for the initial guess.
    //    KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
613

614
615
616

    // === Run Petsc. ===

617
    KSPSolve(solver, petscRhsVec, petscSolVec);
618

619
    // === Transfere values from Petsc's solution vectors to the dof vectors.
620
621
622
623
    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

    for (int i = 0; i < nComponents; i++) {
624
      DOFVector<double> *dofvec = vec.getDOFVector(i);
625
      for (int j = 0; j < nRankDofs; j++)
626
	(*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];      
627
628
629
630
    }

    VecRestoreArray(petscSolVec, &vecPointer);

631
632
633

    // === Synchronize dofs at common dofs, i.e., dofs that correspond to more ===
    // === than one partition.                                                 ===
Thomas Witkowski's avatar
Thomas Witkowski committed
634
    synchVector(vec);
635

636
637
638

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

639
640
641
642
643
644
645
646
647
648
    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);

649
650
651

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

652
    MatDestroy(petscMatrix);
653
654
655
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
    VecDestroy(petscTmpVec);
656
657
658
    KSPDestroy(solver);
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
659
660

  void ParallelDomainBase::synchVector(DOFVector<double> &vec)
661
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
662
    StdMpi<std::vector<double> > stdMpi(mpiComm);
663
664

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
665
	 sendIt != sendDofs.end(); ++sendIt) {
666
      std::vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
667
668
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
669
      
Thomas Witkowski's avatar
Thomas Witkowski committed
670
671
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
672
673
674
675

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

Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
678
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
679

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

Thomas Witkowski's avatar
Thomas Witkowski committed
682
683
684
685
686
    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];
  }
687
688


Thomas Witkowski's avatar
Thomas Witkowski committed
689
690
691
692
693
694
695
696
697
698
699
700
701
702
  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])]);
703
704
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
705
      stdMpi.send(sendIt->first, dofs);
706
707
708
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
712
    stdMpi.startCommunication<double>(MPI_DOUBLE);
713
714

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
715
716
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
717
718

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
719
720
721
722
723
      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++];
724
725
726
727
      }
    }
  }

728
729
730

  void ParallelDomainBase::checkMeshChange()
  {
731
732
    FUNCNAME("ParallelDomainBase::checkMeshChange()");

733
734
735
736
737
    // === 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);
738

739
    if (recvAllValues == 0)
740
741
      return;

742
743
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
744

745
746
747
748
749
    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.
750
751
752
753
754
755
756
757
      RankToBoundMap allBound;
      for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin(); 
	   it != myIntBoundary.boundary.end(); ++it)
	for (std::vector<AtomicBoundary>::iterator bit = it->second.begin();
	     bit != it->second.end(); ++bit)
	  if (bit->rankObj.subObj == FACE)
	    allBound[it->first].push_back(*bit);

758
759
      for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); 
	   it != otherIntBoundary.boundary.end(); ++it)
760
761
762
763
	for (std::vector<AtomicBoundary>::iterator bit = it->second.begin();
	     bit != it->second.end(); ++bit)
	  if (bit->rankObj.subObj == FACE)
	    allBound[it->first].push_back(*bit);
764

765
      // === Check the boundaries and adapt mesh if necessary. ===
766

767
768
769
770
771
772
773
      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);
774
      MSG("LOOP\n");
775
    } while (recvAllValues != 0);
776

777
778

#if (DEBUG != 0)
779
    debug::writeMesh(feSpace, -1, "mesh");
780
781
782
783
784
785
786
787
#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();
788
789
790
791
792
793
794
795
796
797
798
799
800

    // === If there is a non zero pattern computed for Petsc matrix, delete it. So ===
    // === it will be recomputed after next assembling.                            ===

    if (d_nnz) {
      delete [] d_nnz;
      d_nnz = NULL;
    }

    if (o_nnz) {
      delete [] o_nnz;
      o_nnz = NULL;
    }
801
802
803
804
805
806
807
  }

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

808
809
    MSG("CHECK_AND_ADAPT 1!\n");

810
811
812
813
814
    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
815
816
817
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
818
	elCode.init(boundIt->rankObj);
819
820
821
822
	sendCodes[it->first].push_back(elCode);
      }
    }

823
824
    MSG("CHECK_AND_ADAPT 2!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
825
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
826
    stdMpi.send(sendCodes);
827
828
829
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
830
831
832
833
    MSG("CHECK_AND_ADAPT 3!\n");

    clock_t first = clock();

834
    // === Compare received mesh structure codes. ===
835
    
836
    bool meshFitTogether = true;
837
    int superCounter = 0;
838

839
840
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
841
842
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
843
      
844
845
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
846

847
848
849
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
850
851
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
852
853
854
855
856
	  superCounter++;
	  clock_t fff = clock();

	  int refCounter = 0;
	  bool b = fitElementToMeshCode(recvCodes[i], 
857
858
					boundIt->rankObj.el,
					boundIt->rankObj.ithObj, 
859
860
861
862
863
864
865
866
					boundIt->rankObj.elType, refCounter);  

	  MSG("SIZE OF ELINFO3d = %d\n", sizeof(ElInfo3d));
	  MSG("Code length %d with %d refs needed %.5f seconds\n", 
	      recvCodes[i].getNumElements(), 
	      refCounter,
	      TIME_USED(fff, clock()));

867
868
	  if (b)
	    meshFitTogether = false;
869
	}
870
	
871
	i++;
872
873
874
      }
    }

875
876
877
    MSG("time for %d needed %.5f seconds\n", superCounter, TIME_USED(first, clock()));
    MSG("CHECK_AND_ADAPT 4!\n");

878
    return meshFitTogether;
879
  }
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
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


  bool ParallelDomainBase::fitElementToMeshCode(MeshStructure &code, 
						Element *el, 
						int ithSide, 
						int elType,
						int &refCounter)
  {
    FUNCNAME("fitElementToMeshCode()");

    if (code.empty())
      return false;
    
    int s1 = el->getSideOfChild(0, ithSide, elType);
    int s2 = el->getSideOfChild(1, ithSide, elType);

    TEST_EXIT_DBG(s1 != -1 || s2 != -1)("This should not happen!\n");

    if (s1 != -1 && s2 != -1)
      return fitElementToMeshCode2(code, el, ithSide, elType, refCounter);

    if (el->isLeaf()) {
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;
      
      el->setMark(1);
      refineManager->refineElement(el->getMesh(), el);
      refCounter++;
    }
   
    if (s1 != -1)
      return fitElementToMeshCode2(code, 
				   el->getFirstChild(), s1, el->getChildType(elType), refCounter);
    else
      return fitElementToMeshCode2(code, 
				   el->getSecondChild(), s2, el->getChildType(elType), refCounter);
  }

  bool ParallelDomainBase::fitElementToMeshCode2(MeshStructure &code, 
						 Element *el, 
						 int ithSide, 
						 int elType, int &refCounter)
  {
    if (code.isLeafElement())
      return false;

    bool value = false;

    if (el->isLeaf()) {
      el->setMark(1);
      refineManager->refineElement(el->getMesh(), el);     
      value = true;
      refCounter++;
    }

    int s1 = el->getSideOfChild(0, ithSide, elType);
    int s2 = el->getSideOfChild(1, ithSide, elType);
    
    if (s1 != -1) {
      code.nextElement();
      value |= fitElementToMeshCode2(code, el->getFirstChild(), 
				     s1, el->getChildType(elType), refCounter);
    }  

    if (s2 != -1) {
      code.nextElement();	
      value |= fitElementToMeshCode2(code, el->getSecondChild(), 
				     s2, el->getChildType(elType), refCounter);
    }

    return value;
  }

953
954
955
956
  
  void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
  {
    int vecSize = data.size();
957
    SerUtil::serialize(out, vecSize);
958
959
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
960
      SerUtil::serialize(out, dofIndex);
961
962
963
964
965
966
967
968
969
970