Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

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