ParallelDomainBase.cc 92.9 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
#include <boost/lexical_cast.hpp>
5
6
#include "parallel/ParallelDomainBase.h"
#include "parallel/StdMpi.h"
7
8
9
10
11
12
13
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
14
15
#include "DOFMatrix.h"
#include "DOFVector.h"
16
#include "SystemVector.h"
17
#include "VtkWriter.h"
18
#include "ElementDofIterator.h"
19
20
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
21
#include "ElementFileWriter.h"
22
#include "VertexVector.h"
23
#include "MeshStructure.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
24
25
#include "ProblemVec.h"
#include "ProblemInstat.h"
26
#include "Debug.h"
27
28

#include "petscksp.h"
29
30
31

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
32
33
  using boost::lexical_cast;

34
  PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
35
  {    
36
37
    if (iter % 10 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
      std::cout << "[0]  Petsc-Iteration " << iter << ": " << rnorm << std::endl;
38
39
40
41

    return 0;
  }

42
43
44
45
46
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
47
48
49
50
51
52
53
54
55
56
  ParallelDomainBase::ParallelDomainBase(ProblemVec *problemStat,
					 ProblemInstatVec *problemInstat)
    : iterationIF(problemStat),
      timeIF(problemInstat),
      probStat(problemStat),
      name(problemStat->getName()),
      feSpace(problemStat->getFESpace(0)),
      mesh(feSpace->getMesh()),
      refineManager(problemStat->getRefinementManager(0)),
      info(problemStat->getInfo()),
57
      initialPartitionMesh(true),
58
59
      d_nnz(NULL),
      o_nnz(NULL),
60
      nRankDofs(0),
61
      rstart(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
62
      nComponents(problemStat->getNumComponents()),
63
64
      deserialized(false),
      lastMeshChangeIndex(0)
65
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
68
69
70
71
72
    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");

73
74
75
76
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

    // Test if all fe spaces are equal. Yet, different fe spaces are not supported for
    // domain parallelization.
    const FiniteElemSpace *fe = probStat->getFESpace(0);
    for (int i = 0; i < nComponents; i++)
      TEST_EXIT(fe == probStat->getFESpace(i))
	("Parallelization does not supported different FE spaces!\n");

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
    GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
    if (writeSerialization)
      probStat->getFileWriterList().push_back(new Serializer<ParallelDomainBase>(this));

    int readSerialization = 0;
    GET_PARAMETER(0, name + "->input->read serialization", "%d", &readSerialization);
    if (readSerialization) {
      std::string filename = "";
      GET_PARAMETER(0, name + "->input->serialization filename", &filename);
      filename += ".p" + lexical_cast<std::string>(mpiRank);
      MSG("Start serialization with %s\n", filename.c_str());
      std::ifstream in(filename.c_str());
      deserialize(in);
      in.close();
    }
102
103
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
104

105
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
106
  {
107
108
109
110
111
112
113
    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)
114
115
      return;

116
117
118
119
120
    // 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();

121
122
    MSG("START PARTITIONING!\n");

123
124
125
126
127
128
129
    // 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
130
131
#if (DEBUG != 0)
    ElementIdxToDofs elMap;
132
    dbgCreateElementMap(elMap);
133
134
135
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
136

137
138
139
140
141
      if (writePartMesh > 0)
	writePartitioningMesh("part.vtu");
      else 
	MSG("Skip write part mesh!\n");
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
142
#endif
143

144
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
145

146
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
147

148
149
150
151
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
152
153
154
155
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

156
157
158
159
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

160
161
    // === If in debug mode, make some tests. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
162
#if (DEBUG != 0)
163
    MSG("AMDiS runs in debug mode, so make some test ...\n");
164
165
    dbgTestElementMap(elMap);
    dbgTestInteriorBoundary();
166
    dbgTestCommonDofs(true);
167
    MSG("Debug mode tests finished!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
168

169
170
    debug::writeMesh(feSpace, -1, "macromesh");   
#endif
171

172
173
174
    // === Create periodic dof mapping, if there are periodic boundaries. ===

    createPeriodicMap();
175

176
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
177

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

Thomas Witkowski's avatar
Thomas Witkowski committed
181
    if (globalRefinement > 0) {
182
      refineManager->globalRefine(mesh, globalRefinement);
183

184
      updateLocalGlobalNumbering();
185

186
      // === Update periodic mapping, if there are periodic boundaries. ===
187

188
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
189
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229


    /// === Set DOF rank information to all matrices and vectors. ===

    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++)
 	if (probStat->getSystemMatrix(i, j))
 	  probStat->getSystemMatrix(i, j)->setRankDofs(isRankDof);      

      TEST_EXIT_DBG(probStat->getRHS()->getDOFVector(i))("No rhs vector!\n");
      TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i))("No solution vector!\n");

      probStat->getRHS()->getDOFVector(i)->setRankDofs(isRankDof);
      probStat->getSolution()->getDOFVector(i)->setRankDofs(isRankDof);
    }

    // === Remove periodic boundary conditions in sequential problem definition. ===

    // Remove periodic boundaries in boundary manager on matrices and vectors.
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	DOFMatrix* mat = probStat->getSystemMatrix(i, j);
 	if (mat && mat->getBoundaryManager())
	  removeBoundaryCondition(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
      }

      if (probStat->getSolution()->getDOFVector(i)->getBoundaryManager())
	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat->getSolution()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));

      if (probStat->getRHS()->getDOFVector(i)->getBoundaryManager())
	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat->getRHS()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));
    }

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh,  -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
230
231
  }

232
233

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
234
  {}
235

236
237
238
  
  void ParallelDomainBase::updateDofAdmins()
  {
239
240
241
    FUNCNAME("ParallelDomainBase::updateDofAdmins()");

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
242
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
243
244
245

      // There must be always more allocated DOFs than used DOFs in DOFAdmin. Otherwise,
      // it is not possible to define a first DOF hole.
246
      if (static_cast<unsigned int>(admin.getSize()) == mapLocalGlobalDofs.size())
247
	admin.enlargeDOFLists();
248
249
250
      
      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
251
      for (unsigned int j = 0; j < mapLocalGlobalDofs.size(); j++)
252
253
 	admin.setDOFFree(j, false);

254
255
256
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
257
258
259
    }
  }

260

261
262
263
264
265
266
267
268
269
270
  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)
271
	("Mesh is already refined! This does not work with parallelization!\n");
272
273
274
275
276
277
278
279
280
281
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

282

283
284
  void ParallelDomainBase::setDofMatrix(DOFMatrix* mat, int dispMult, 
					int dispAddRow, int dispAddCol)
285
  {
286
287
288
289
    FUNCNAME("ParallelDomainBase::setDofMatrix()");

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

290
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
291
292
293
294
295
296
    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());

297
    typedef traits::range_generator<row, Matrix>::type cursor_type;
298
299
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

300
301
302
303
304
    std::vector<int> cols;
    std::vector<double> values;
    cols.reserve(300);
    values.reserve(300);

305
306
307
    // === Traverse all rows of the dof matrix and insert row wise the values ===
    // === to the petsc matrix.                                               ===

308
309
310
311
312
313
    for (cursor_type cursor = begin<row>(mat->getBaseMatrix()), 
	   cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {

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

314
      // Global index of the current row dof.
315
      int globalRowDof = mapLocalGlobalDofs[*cursor];
316
      // Test if the current row dof is a periodic dof.
317
      bool periodicRow = (periodicDof.count(globalRowDof) > 0);
318

319
320
321

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

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

	// Set only non null values.
328
	if (value(*icursor) != 0.0) {
329
	  // Global index of the current column index.
330
	  int globalColDof = mapLocalGlobalDofs[col(*icursor)];
331
	  // Calculate the exact position of the column index in the petsc matrix.
332
333
	  int colIndex = globalColDof * dispMult + dispAddCol;

334
335
	  // 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.
336
 	  if (periodicRow == false && periodicDof.count(globalColDof) > 0) {
337
338
339
340
341
	    // 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.
342
 	    cols.push_back(colIndex);
343
 	    values.push_back(value(*icursor) * scalFactor);
344

345
346
347
	    // Insert the periodic entries.
 	    for (std::set<DegreeOfFreedom>::iterator it = 
		   periodicDof[globalColDof].begin();
348
349
 		 it != periodicDof[globalColDof].end(); ++it) {
 	      cols.push_back(*it * dispMult + dispAddCol);
350
351
 	      values.push_back(value(*icursor) * scalFactor);
	    }
352
 	  } else {
353
	    // Neigher row nor column dof index is periodic, simple add entry.
354
355
356
	    cols.push_back(colIndex);
	    values.push_back(value(*icursor));
	  }
357
	}
358
359
      }

360
361
362
363
364

      // === 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.
365
      int rowIndex = globalRowDof * dispMult + dispAddRow;
366
      
367
      if (periodicRow) {
368
369
370
	// The row dof is periodic, so send dof to all the corresponding rows.

	double scalFactor = 1.0 / (periodicDof[globalRowDof].size() + 1.0);
371
	
372
	int diagIndex = -1;
373
	for (int i = 0; i < static_cast<int>(values.size()); i++) {
374
375
376
377
378
	  // 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;
379
380
381
	  else
	    diagIndex = i;
	}
382
383
384
385
386
387
388
	
	// 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.
389
390
391
	if (diagIndex != -1)
	  values[diagIndex] = 0.0;

392
	// Send the row to all periodic row indices.
393
394
395
	for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRowDof].begin();
	     it != periodicDof[globalRowDof].end(); ++it) {
	  int perRowIndex = *it * dispMult + dispAddRow;
396
397
	  MatSetValues(petscMatrix, 1, &perRowIndex, cols.size(), 
		       &(cols[0]), &(values[0]), ADD_VALUES);
398
399
400
	}

      } else {
401
402
403
	// 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);
404
      }    
405
    }
406
  }
407

408

409
  void ParallelDomainBase::setDofVector(Vec& petscVec, DOFVector<double>* vec, 
410
411
					int dispMult, int dispAdd)
  {
412
    // Traverse all used dofs in the dof vector.
413
414
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
415
      // Calculate global row index of the dof.
416
      int globalRow = mapLocalGlobalDofs[dofIt.getDOFIndex()];
417
      // Calculate petsc index of the row dof.
418
419
420
      int index = globalRow * dispMult + dispAdd;

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

	double value = *dofIt / (periodicDof[globalRow].size() + 1.0);
424
425
426
427
428
429
430
	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);
	}
431

432
      } else {
433
	// The dof index is not periodic.
434
435
436
	double value = *dofIt;
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
      }
437
    }    
438
439
  }

440

441
  void ParallelDomainBase::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
442
  {
443
    FUNCNAME("ParallelDomainBase::createPetscNnzStructure()");
444

445
446
    TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n");
    TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n");
447

448
449
450
451
452
453
    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;
    }
454

455
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
456
    namespace traits = mtl::traits;
457
    typedef DOFMatrix::base_matrix_type Matrix;
458
    typedef std::vector<std::pair<int, int> > MatrixNnzEntry;
459

460
461
462
463
464
465
    // 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;

466
467
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
468
469
470
471
472
473
 	if ((*mat)[i][j]) {
	  Matrix bmat = (*mat)[i][j]->getBaseMatrix();

	  traits::col<Matrix>::type col(bmat);
	  traits::const_value<Matrix>::type value(bmat);
	  
474
	  typedef traits::range_generator<row, Matrix>::type cursor_type;
475
476
	  typedef traits::range_generator<nz, cursor_type>::type icursor_type;
	  
477
478
479
	  for (cursor_type cursor = begin<row>(bmat), 
		 cend = end<row>(bmat); cursor != cend; ++cursor) {

480
481
482
	    // 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;
483

484
	    if (isRankDof[*cursor]) {
485
486
487
488
489
490

	      // === 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;
491
492

#if (DEBUG != 0)    
493
	      if (localPetscRowIdx < 0 || localPetscRowIdx >= nRankRows) {
494
		std::cout << "ERROR in rank: " << mpiRank << std::endl;
495
496
		std::cout << "  Wrong r = " << localPetscRowIdx << " " << *cursor 
			  << " " << mapLocalGlobalDofs[*cursor] << " " 
497
498
499
500
			  << nRankRows << std::endl;
		ERROR_EXIT("Should not happen!\n");
	      }
#endif
501
	      
502
	      // Traverse all non zero entries in this row.
503
504
505
	      for (icursor_type icursor = begin<nz>(cursor), 
		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
		if (value(*icursor) != 0.0) {
506
		  int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
507

508
509
510
511
512
		  // 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]++;
513
		  else
514
		    o_nnz[localPetscRowIdx]++;
515
516
517
		}    
	      }
	    } else {
518
519
520
521
522
	      
	      // === 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.                          ===
523

524
525
	      // Find out who is the member of this DOF.
	      int sendToRank = -1;
526
527
528
529
530
	      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) {
531
532
533
534
535
536

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

537
538
539
540
541
542
543
544
545
546
547
		    sendToRank = it->first;
		    break;
		  }
		}

		if (sendToRank != -1)
		  break;
	      }

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

548
	      // Send all non zero entries to the member of the row DOF.
549
550
551
	      for (icursor_type icursor = begin<nz>(cursor), 
		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
		if (value(*icursor) != 0.0) {
552
		  int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
553
		  
554
555
		  sendMatrixEntry[sendToRank].
		    push_back(std::make_pair(petscRowIdx, petscColIdx));
556
557
558
		}
	      }

559
560
561
562
	    } // if (isRankDof[*cursor]) ... else ...
	  } // for each row in mat[i][j]
	} // if mat[i][j]
      } 
563
564
    }

565
    // === Send and recv the nnz row structure to/from other ranks. ===
566

Thomas Witkowski's avatar
Thomas Witkowski committed
567
    StdMpi<MatrixNnzEntry> stdMpi(mpiComm, true);
568
569
570
    stdMpi.send(sendMatrixEntry);
    stdMpi.recv(sendDofs);
    stdMpi.startCommunication<int>(MPI_INT);
Thomas Witkowski's avatar
Thomas Witkowski committed
571

572
573
    // === Evaluate the nnz structure this rank got from other ranks and add it to ===
    // === the PETSc nnz data structure.                                           ===
574

575
576
577
578
579
580
    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;
581

582
	  int localRowIdx = r - rstart * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
583

584
585
586
	  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);
587
	  
588
	  if (c < rstart * nComponents || c >= rstart * nComponents + nRankRows)
589
	    o_nnz[localRowIdx]++;
590
	  else
591
	    d_nnz[localRowIdx]++;
592
593
	}
      }
594
    }  
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
  }


  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
620

621
    // === Create PETSc matrix with the computed nnz data structure. ===
622
623

    MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
624
625
		    0, d_nnz, 0, o_nnz, &petscMatrix);

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

628
629
630
631
632
633
#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
634

635
    // === Transfer values from DOF matrices to the PETSc matrix. === 
636
637
638

    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
639
	if ((*mat)[i][j])
640
641
	  setDofMatrix((*mat)[i][j], nComponents, i, j);

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

644
645
646
    MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);

647
648
    // === Transfer values from DOF vector to the PETSc vector. === 

649
    for (int i = 0; i < nComponents; i++)
650
      setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i);
651

652
653
654
    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);

655
    INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", TIME_USED(first, clock()));
656
657
658
  }


659
  void ParallelDomainBase::solvePetscMatrix(SystemVector &vec)
660
661
662
  {
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

663
664
665
666
667
668
669
670
671
#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

672
    // === Init Petsc solver. ===
673

674
    KSP solver;
675
676
    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 
677
    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
678
679
    KSPSetType(solver, KSPBCGS);
    KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0);
680
    KSPSetFromOptions(solver);
681
682
    // Do not delete the solution vector, use it for the initial guess.
    //    KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
683

684
685
686

    // === Run Petsc. ===

687
    KSPSolve(solver, petscRhsVec, petscSolVec);
688

689
    // === Transfere values from Petsc's solution vectors to the dof vectors.
690
691
692
693
    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

    for (int i = 0; i < nComponents; i++) {
694
      DOFVector<double> *dofvec = vec.getDOFVector(i);
695
      for (int j = 0; j < nRankDofs; j++)
696
	(*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];      
697
698
699
700
    }

    VecRestoreArray(petscSolVec, &vecPointer);

701
702
703

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

706
707
708

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

709
710
711
712
713
714
715
716
717
718
    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);

719
720
721

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

722
    MatDestroy(petscMatrix);
723
724
725
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
    VecDestroy(petscTmpVec);
726
727
728
    KSPDestroy(solver);
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
729
730

  void ParallelDomainBase::synchVector(DOFVector<double> &vec)
731
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
732
    StdMpi<std::vector<double> > stdMpi(mpiComm);
733
734

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
735
	 sendIt != sendDofs.end(); ++sendIt) {
736
      std::vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
737
738
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
739
      
Thomas Witkowski's avatar
Thomas Witkowski committed
740
741
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
742
743
744
745

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

Thomas Witkowski's avatar
Thomas Witkowski committed
746
747
748
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
749

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

Thomas Witkowski's avatar
Thomas Witkowski committed
752
753
754
755
756
    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];
  }
757
758


Thomas Witkowski's avatar
Thomas Witkowski committed
759
760
761
762
763
764
765
766
767
768
769
770
771
772
  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])]);
773
774
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
775
      stdMpi.send(sendIt->first, dofs);
776
777
778
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
782
    stdMpi.startCommunication<double>(MPI_DOUBLE);
783
784

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
785
786
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
787
788

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
789
790
791
792
793
      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++];
794
795
796
797
      }
    }
  }

798

Thomas Witkowski's avatar
Thomas Witkowski committed
799
800
801
802
803
804
805
806
807
808
809
810
  void ParallelDomainBase::removeBoundaryCondition(BoundaryIndexMap& boundaryMap)
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


811
812
  void ParallelDomainBase::checkMeshChange()
  {
813
814
    FUNCNAME("ParallelDomainBase::checkMeshChange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
815
816
    int dim = mesh->getDim();

817
818
819
820
821
    // === 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);
822

823
    if (recvAllValues == 0)
824
825
      return;

826
827
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
828

829
830
831
832
833
    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.
834
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
835
836
837
838
839
840
841
842
843

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
	if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
	  allBound[it.getRank()].push_back(it.getBound());

      for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
	if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
	  allBound[it.getRank()].push_back(it.getBound());

844

845
      // === Check the boundaries and adapt mesh if necessary. ===
846

847
848
849
850
851
852
853
854
      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);
855

856
857

#if (DEBUG != 0)
858
    debug::writeMesh(feSpace, -1, "mesh");
859
860
861
862
863
864
865
866
#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();
867
868
869
870
871
872
873
874
875
876
877
878
879

    // === 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;
    }
880
881
882
883
884
885
886
887
888
889
890
891
  }

  
  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) {
892
893
894
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
895
	elCode.init(boundIt->rankObj);
896
897
898
899
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
900
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
901
    stdMpi.send(sendCodes);
902
903
904
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
905
    // === Compare received mesh structure codes. ===
906
    
907
908
    bool meshFitTogether = true;

909
910
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
911
912
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
913
      
914
915
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
916

917
918
919
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
920
921
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
922
923
924

	  int refCounter = 0;
	  bool b = fitElementToMeshCode(recvCodes[i], 
925
926
					boundIt->rankObj.el,
					boundIt->rankObj.ithObj, 
927
928
					boundIt->rankObj.elType, refCounter);  

929
930
	  if (b)
	    meshFitTogether = false;
931
	}
932
	
933
	i++;
934
935
936
      }
    }

937
    return meshFitTogether;
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


  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