ParallelDebug.cc 28.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
#include "parallel/ParallelDebug.h"
#include "parallel/MeshDistributor.h"
#include "parallel/MpiHelper.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
16
17
18
19
20
#include "ProblemVec.h"
#include "DOFVector.h"
#include "FixVec.h"
#include "StdMpi.h"
#include "Debug.h"
21
#include "io/VtkWriter.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
22
23
24

namespace AMDiS {

25
  void ParallelDebug::testInteriorBoundary(MeshDistributor &pdb)
Thomas Witkowski's avatar
Thomas Witkowski committed
26
  {
27
    FUNCNAME("ParallelDebug::testInteriorBoundary()");
Thomas Witkowski's avatar
Thomas Witkowski committed
28

29
    typedef MeshDistributor::RankToBoundMap RankToBoundMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
30
31
32
33

    std::vector<int*> sendBuffers, recvBuffers;

    MPI::Request request[pdb.myIntBoundary.boundary.size() + 
34
35
			 pdb.otherIntBoundary.boundary.size() +
                         pdb.periodicBoundary.boundary.size() * 2];
Thomas Witkowski's avatar
Thomas Witkowski committed
36
37
    int requestCounter = 0;

38
39
40

    // === Send rank's boundary information. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
    for (RankToBoundMap::iterator rankIt = pdb.myIntBoundary.boundary.begin();
	 rankIt != pdb.myIntBoundary.boundary.end(); ++rankIt) {

      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
      for (int i = 0; i < nSendInt; i++)
	buffer[i] = (rankIt->second)[i].rankObj.elIndex;
      
      sendBuffers.push_back(buffer);
      
      request[requestCounter++] =
	pdb.mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
    }

55
56
57

    // === Receive information from other ranks about the interior boundaries. ====

Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
60
61
62
63
64
65
66
67
68
    for (RankToBoundMap::iterator rankIt = pdb.otherIntBoundary.boundary.begin();
	 rankIt != pdb.otherIntBoundary.boundary.end(); ++rankIt) {
      int nRecvInt = rankIt->second.size();
      int *buffer = new int[nRecvInt];
      recvBuffers.push_back(buffer);

      request[requestCounter++] = 
	pdb.mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
    }


69
70
    // === To the last, do the same of periodic boundaries. ===

71
    for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();       
72
	 rankIt != pdb.periodicBoundary.boundary.end(); ++rankIt) {
73
74
75
      if (rankIt->first == pdb.mpiRank)
	continue;

76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
      int nValues = rankIt->second.size();
      int* sBuffer = new int[nValues];
      for (int i = 0; i < nValues; i++)
	sBuffer[i] = (rankIt->second)[i].rankObj.elIndex;

      sendBuffers.push_back(sBuffer);

      request[requestCounter++] =
	pdb.mpiComm.Isend(sBuffer, nValues, MPI_INT, rankIt->first, 0);

      int *rBuffer = new int[nValues];
      recvBuffers.push_back(rBuffer);

      request[requestCounter++] = 
	pdb.mpiComm.Irecv(rBuffer, nValues, MPI_INT, rankIt->first, 0);      
    }

    // === Finish communication and delete all send buffers. ===

    MPI::Request::Waitall(requestCounter, request);
Thomas Witkowski's avatar
Thomas Witkowski committed
96
97
98
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i];

99
100
101
102
103

    // === Finally, check the results, i.e., the indices of element at the     === 
    // === boundaries, if they fit together. First check the interior bounds,  ===
    // === and after this the periodic ones.                                   ===

Thomas Witkowski's avatar
Thomas Witkowski committed
104
105
106
107
    int bufCounter = 0;
    for (RankToBoundMap::iterator rankIt = pdb.otherIntBoundary.boundary.begin();
	 rankIt != pdb.otherIntBoundary.boundary.end(); ++rankIt) {

108
109
      TEST_EXIT(rankIt->second.size() == 
		pdb.otherIntBoundary.boundary[rankIt->first].size())
Thomas Witkowski's avatar
Thomas Witkowski committed
110
111
	("Boundaries does not fit together!\n");      

112
      for (unsigned int i = 0; i < rankIt->second.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
113
114
115
116
117
118
119
120
	int elIndex1 = recvBuffers[bufCounter][i];
	int elIndex2 = pdb.otherIntBoundary.boundary[rankIt->first][i].neighObj.elIndex;

	TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at interior boundary!\n");
      }

      delete [] recvBuffers[bufCounter++];
    }
121
122
123
124


    for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();
	 rankIt != pdb.periodicBoundary.boundary.end(); ++rankIt) {
125
126
127
      if (rankIt->first == pdb.mpiRank)
	continue;

128
129
130
131
      for (unsigned int i = 0; i < rankIt->second.size(); i++) {
	int elIndex1 = recvBuffers[bufCounter][i];
	int elIndex2 = pdb.periodicBoundary.boundary[rankIt->first][i].neighObj.elIndex;

132
133
134
135
	TEST_EXIT(elIndex1 == elIndex2)
	  ("Wrong element index at periodic boundary el %d with rank %d: %d %d\n", 
	   pdb.periodicBoundary.boundary[rankIt->first][i].rankObj.elIndex,
	   rankIt->first, elIndex1, elIndex2);
136
137
138
139
      }

      delete [] recvBuffers[bufCounter++];
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
140
141
142
  }


143
144
145
146
  void ParallelDebug::testPeriodicBoundary(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testPeriodicBoundary()");

147
148
    return;

149
150
    // === 1. check: All periodic DOFs should have at least a correct number ===
    // === of periodic associations.                                         ===
151
   
152
153
154
155
156
157
158
159
160
161
162
163
164
165
    for (map<int, std::set<BoundaryType> >::iterator it = pdb.periodicDofAssociations.begin();
	 it != pdb.periodicDofAssociations.end(); ++it) {
      WorldVector<double> c;
      pdb.mesh->getDofIndexCoords(it->first, pdb.feSpace, c);
      int nAssoc = it->second.size();
      TEST_EXIT_DBG(nAssoc == 1 || nAssoc == 3 || (pdb.mesh->getDim() == 3 && nAssoc == 7))
	("Should not happen! DOF %d (%e %e %e) has %d periodic associations!\n", 
	 it->first, c[0], c[1], (pdb.mesh->getDim() == 2 ? 0.0 : c[2]), nAssoc);
    }    


    // === 2. check: All periodic DOFs must be symmetric, i.e., if A is mapped ===
    // === to B, then B must be mapped to A.                                   ===

166
167
168
169
170
171
172
173
174
175
176
    typedef MeshDistributor::PeriodicDofMap PeriodicDofMap;

    StdMpi<PeriodicDofMap> stdMpi(pdb.mpiComm, true);

    if (pdb.mpiRank == 0) {
      for (int i = 1; i < pdb.mpiSize; i++)
	stdMpi.recv(i);
    } else {
      stdMpi.send(0, pdb.periodicDof);
    }

177
    stdMpi.startCommunication();
178
179
180
181
182
183

    int foundError = 0;

    // === The boundary DOFs are checked only on the zero rank. === 

    if (pdb.mpiRank == 0) {
184
      // Stores to each rank the periodic DOF mappings of this rank.
185
186
187
188
189
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
230
231
232
233
234
235
236
237
238
239
240
241
242
      std::map<int, PeriodicDofMap> rankToMaps;
      PeriodicDofMap dofMap = pdb.periodicDof;
      rankToMaps[0] = dofMap;

      for (int i = 1; i < pdb.mpiSize; i++) {
	PeriodicDofMap &otherMap = stdMpi.getRecvData(i);
	rankToMaps[i] = otherMap;
	
	for (PeriodicDofMap::iterator it = otherMap.begin(); 
	     it != otherMap.end(); ++it) {
	  for (MeshDistributor::DofMapping::iterator dofIt = it->second.begin();
	       dofIt != it->second.end(); ++dofIt) {
	    if (dofMap.count(it->first) == 1 &&
		dofMap[it->first].count(dofIt->first) == 1) {
	      TEST_EXIT_DBG(dofMap[it->first][dofIt->first] == dofIt->second)
		("Should not happen!\n");
	    } else {
	      dofMap[it->first][dofIt->first] = dofIt->second;
	    }
	  }
	}
      }


      // === Now we test if global DOF A is mapped to B, then B must be mapped ===
      // === to A for the same boundary type.                                  ===

      for (PeriodicDofMap::iterator it = dofMap.begin(); 
	   it != dofMap.end(); ++it) {
	for (MeshDistributor::DofMapping::iterator dofIt = it->second.begin();
	     dofIt != it->second.end(); ++dofIt) {
	  if (it->second[dofIt->second] != dofIt->first) {
	    MSG("[DBG]  For boundary type %d: DOF %d -> %d, but %d -> %d!\n",
		it ->first, 
		dofIt->first, dofIt->second, 
		dofIt->second, it->second[dofIt->second]);

	    for (int i = 0; i < pdb.mpiSize; i++) {
	      if (rankToMaps[i][it->first].count(dofIt->first) == 1) {
		MSG("[DBG]    %d -> %d in rank %d\n", 
		    dofIt->first, rankToMaps[i][it->first][dofIt->first], i);
	      }

	      if (rankToMaps[i][it->first].count(dofIt->second) == 1) {
		MSG("[DBG]    %d -> %d in rank %d\n", 
		    dofIt->second, rankToMaps[i][it->first][dofIt->second], i);
	      }
	    }
	    
	    ERROR("Wrong periodic DOFs!\n");
	    foundError = 1;
	  }
	}
      }
    }

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262


    // === 3. check: On all edge and face periodic DOFs, at least on coordinate of ===
    // === each periodic DOF pair must be equal (at least as long we consider      ===
    // === periodic boundaries only on rectangulars and boxes.                     ===

    RankToCoords sendCoords;
    std::map<int, std::vector<BoundaryType> > rankToDofType;

    for (InteriorBoundary::RankToBoundMap::iterator it = pdb.periodicBoundary.boundary.begin();
	 it != pdb.periodicBoundary.boundary.end(); ++it) {
      if (it->first == pdb.mpiRank)
	continue;

      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {	
	if (boundIt->rankObj.subObj == VERTEX)
	  continue;

	DofContainer dofs;
263
	boundIt->rankObj.el->getAllDofs(pdb.feSpace, boundIt->rankObj, dofs);
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
	
	for (unsigned int i = 0; i < dofs.size(); i++) {
	  WorldVector<double> c;
	  pdb.mesh->getDofIndexCoords(*(dofs[i]), pdb.feSpace, c);
	  sendCoords[it->first].push_back(c);
	  rankToDofType[it->first].push_back(boundIt->type);
	}
      }
    }

    // Each rank must receive exactly the same number of coordinates as it sends
    // to another rank.
    RankToCoords recvCoords;
    for (RankToCoords::iterator it = sendCoords.begin(); 
	 it != sendCoords.end(); ++it)
      recvCoords[it->first].resize(it->second.size());


    StdMpi<CoordsVec> stdMpiCoords(pdb.mpiComm, true);
    stdMpiCoords.send(sendCoords);
    stdMpiCoords.recv(recvCoords);   
285
    stdMpiCoords.startCommunication();
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315


    for (RankToCoords::iterator it = sendCoords.begin(); 
	 it != sendCoords.end(); ++it) {
      for (unsigned int i = 0; i < it->second.size(); i++) {
	WorldVector<double> &c0 = it->second[i];
	WorldVector<double> &c1 = stdMpiCoords.getRecvData(it->first)[i];


	int nEqual = 0;
	for (int j = 0; j < pdb.mesh->getDim(); j++)
	  if (c0[j] == c1[j])
	    nEqual++;

	if ((rankToDofType[it->first][i] >= -3 && nEqual < 2) ||
	    (rankToDofType[it->first][i] < -3 && nEqual == 0)) {
	  MSG("[DBG]  %d-ith periodic DOF in boundary between ranks %d <-> %d is not correct!\n",
	      i, pdb.mpiRank, it->first);
	  MSG("[DBG]  Coords on rank %d: %f %f %f\n", 
	      pdb.mpiRank, c0[0], c0[1], (pdb.mesh->getDim() == 3 ? c0[2] : 0.0));
	  MSG("[DBG]  Coords on rank %d: %f %f %f\n", 
	      it->first, c1[0], c1[1], (pdb.mesh->getDim() == 3 ? c1[2] : 0.0));

	  foundError = 1;
	}
      }
    }

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Wrond periodic coordinates found on at least on rank!\n");
316
317
318
  }


319
  void ParallelDebug::testCommonDofs(MeshDistributor &pdb, bool printCoords)
Thomas Witkowski's avatar
Thomas Witkowski committed
320
  {
321
    FUNCNAME("ParallelDebug::testCommonDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
324
325

    clock_t first = clock();

    int testCommonDofs = 1;
326
    Parameters::get("dbg->test common dofs", testCommonDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
327
328
329
330
331
332
333
334
335
    if (testCommonDofs == 0) {
      MSG("Skip test common dofs!\n");
      return;
    }

    /// Defines a mapping type from rank numbers to sets of DOFs.
    typedef std::map<int, DofContainer> RankToDofContainer;

    // Maps to each neighbour rank an array of WorldVectors. This array contains the 
336
    // coordinates of all DOFs this rank shares on the interior boundary with the 
Thomas Witkowski's avatar
Thomas Witkowski committed
337
    // neighbour rank. A rank sends the coordinates to another rank, if it owns the
338
    // boundarys DOFs.
Thomas Witkowski's avatar
Thomas Witkowski committed
339
340
    RankToCoords sendCoords;

341
342
    // A rank receives all boundary DOFs that are at its interior boundaries but are
    // not owned by the rank. This map stores for each rank the coordinates of DOFs
Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
345
346
347
    // this rank expectes to receive from.
    RankToCoords recvCoords;

    DOFVector<WorldVector<double> > coords(pdb.feSpace, "dofCorrds");
    pdb.mesh->getDofIndexCoords(pdb.feSpace, coords);
348

Thomas Witkowski's avatar
Thomas Witkowski committed
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
    for (RankToDofContainer::iterator it = pdb.sendDofs.begin();
	 it != pdb.sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin(); 
	   dofIt != it->second.end(); ++dofIt)
	sendCoords[it->first].push_back(coords[**dofIt]);

    for (RankToDofContainer::iterator it = pdb.recvDofs.begin();
	 it != pdb.recvDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	recvCoords[it->first].push_back(coords[**dofIt]);

    std::vector<int> sendSize(pdb.mpiSize, 0);
    std::vector<int> recvSize(pdb.mpiSize, 0);
    std::vector<int> recvSizeBuffer(pdb.mpiSize, 0);
    MPI::Request request[(pdb.mpiSize - 1) * 2];
    int requestCounter = 0;

    for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it)
      sendSize[it->first] = it->second.size();

    for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it)
      recvSize[it->first] = it->second.size();

    for (int i = 0; i < pdb.mpiSize; i++) {
      if (i == pdb.mpiRank)
	continue;

377
378
      request[requestCounter++] = 
	pdb.mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
381
382
383
384
    }   

    for (int i = 0; i < pdb.mpiSize; i++) {
      if (i == pdb.mpiRank)
	continue;

385
386
      request[requestCounter++] = 
	pdb.mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
387
388
389
390
391
    }

    MPI::Request::Waitall(requestCounter, request);


392
    int foundError = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
393
394
395
396
397
    for (int i = 0; i < pdb.mpiSize; i++) {
      if (i == pdb.mpiRank)
	continue;

      if (recvSize[i] != recvSizeBuffer[i]) {
398
399
400
	ERROR("MPI rank %d expectes to receive %d DOFs from rank %d. But this rank sends %d DOFs!\n", 
	      pdb.mpiRank, recvSize[i], i, recvSizeBuffer[i]);	
	foundError = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
401
402
      }
    }
403
404
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
407
408
409
410
411
412

    // === Now we know that the number of send and received DOFs fits together. ===
    // === So we can check if also the coordinates of the communicated DOFs are ===
    // === the same on both corresponding ranks.                                ===

    StdMpi<CoordsVec> stdMpi(pdb.mpiComm, true);
    stdMpi.send(sendCoords);
    stdMpi.recv(recvCoords);   
413
    stdMpi.startCommunication();
Thomas Witkowski's avatar
Thomas Witkowski committed
414
415
416
417
418
419
420

    int dimOfWorld = Global::getGeo(WORLD);

    // === Compare the received with the expected coordinates. ===

    for (RankToCoords::iterator it = stdMpi.getRecvData().begin(); 
	 it != stdMpi.getRecvData().end(); ++it) {
421
422
423
      for (unsigned int i = 0; i < it->second.size(); i++) {
	WorldVector<double> tmp = (it->second)[i];
	tmp -=  recvCoords[it->first][i];
424

425
426
427
428
	if (norm(tmp) > 1e-13) {
	  // === Print error message if the coordinates are not the same. ===
	  if (printCoords) {
	    MSG("[DBG] i = %d\n", i);	  
429
430
431
	    std::stringstream oss;
	    oss.precision(5);
	    oss << "[DBG] Rank " << pdb.mpiRank << " from rank " << it->first
432
		<< " expect coords (";
Thomas Witkowski's avatar
Thomas Witkowski committed
433
	    for (int k = 0; k < dimOfWorld; k++) {
434
	      oss << recvCoords[it->first][i][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
435
	      if (k + 1 < dimOfWorld)
436
		oss << " / ";
Thomas Witkowski's avatar
Thomas Witkowski committed
437
	    }
438
	    oss << ")  received coords (";
Thomas Witkowski's avatar
Thomas Witkowski committed
439
	    for (int k = 0; k < dimOfWorld; k++) {
440
	      oss << (it->second)[i][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
441
	      if (k + 1 < dimOfWorld)
442
		oss << " / ";
Thomas Witkowski's avatar
Thomas Witkowski committed
443
	    }
444
445
	    oss << ")";
	    MSG("%s\n", oss.str().c_str());
446
	    
Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
	    debug::printInfoByDof(pdb.feSpace, *(pdb.recvDofs[it->first][i]));
	  }
449
450
451
452
	  ERROR("Wrong DOFs in rank %d!\n", pdb.mpiRank);
	  foundError = 1;
	}	 
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
453
    }
454
455
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
458
459
460

    INFO(pdb.info, 8)("Test common dofs needed %.5f seconds\n", TIME_USED(first, clock()));
  }


461
462
463
464
465
466
467
468
469
470
471
472
  void ParallelDebug::testGlobalIndexByCoords(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testGlobalIndexByCoords()");

    DOFVector<WorldVector<double> > coords(pdb.feSpace, "tmp");
    pdb.mesh->getDofIndexCoords(pdb.feSpace, coords);

    typedef std::map<WorldVector<double>, int> CoordsIndexMap;
    CoordsIndexMap coordsToIndex;

    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it)
473
      coordsToIndex[(*it)] = pdb.mapLocalGlobalDofs[it.getDOFIndex()];    
474
475
476
477
478
479
480
481
482

    StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true);
    for (RankToDofContainer::iterator it = pdb.sendDofs.begin();
	 it != pdb.sendDofs.end(); ++it)
      stdMpi.send(it->first, coordsToIndex);
    for (RankToDofContainer::iterator it = pdb.recvDofs.begin();
	 it != pdb.recvDofs.end(); ++it)
      stdMpi.recv(it->first);
   
483
    stdMpi.startCommunication();
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513

    int foundError = 0;
    for (RankToDofContainer::iterator it = pdb.recvDofs.begin();
	 it != pdb.recvDofs.end(); ++it) {
      CoordsIndexMap& otherCoords = stdMpi.getRecvData(it->first);

      for (CoordsIndexMap::iterator coordsIt = otherCoords.begin();
	   coordsIt != otherCoords.end(); ++coordsIt) {
	if (coordsToIndex.count(coordsIt->first) == 1 &&
	    coordsToIndex[coordsIt->first] != coordsIt->second) {
	  std::stringstream oss;
	  oss.precision(5);
	  oss << "DOF at coords ";
	  for (int i = 0; i < Global::getGeo(WORLD); i++)
	    oss << coordsIt->first[i] << " ";
	  oss << " do not fit together on rank " 
	      << pdb.getMpiRank() << " (global index: " 
	      << coordsToIndex[coordsIt->first] << " and on rank "
	      << it->first << " (global index: " << coordsIt->second << ")";

	  MSG("[DBG] %s\n", oss.str().c_str());
	  foundError = 1;
	}
      }
    }

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
  }

514

515
  void ParallelDebug::testAllElements(MeshDistributor &pdb)
516
  {
517
    FUNCNAME("ParallelDebug::testAllElements()");
518

519
520
521
522
    std::set<int> macroElements;
    int minElementIndex = std::numeric_limits<int>::max();
    int maxElementIndex = std::numeric_limits<int>::min();   

523
    TraverseStack stack;
524
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, 0, Mesh::CALL_EL_LEVEL);
525
    while (elInfo) {
526
527
528
529
530
531
532
533
534
535
      int elIndex = elInfo->getElement()->getIndex();
      minElementIndex = std::min(minElementIndex, elIndex);
      maxElementIndex = std::max(maxElementIndex, elIndex);
      macroElements.insert(elInfo->getElement()->getIndex());
      elInfo = stack.traverseNext(elInfo);      
    }
    
    int globalMinIndex, globalMaxIndex;
    pdb.mpiComm.Allreduce(&minElementIndex, &globalMinIndex, 1, MPI_INT, MPI_MIN);
    pdb.mpiComm.Allreduce(&maxElementIndex, &globalMaxIndex, 1, MPI_INT, MPI_MAX);
536

537
538
539
    TEST_EXIT(globalMinIndex == 0)("No macro element with index 0!\n");
    for (int i = 0; i <= globalMaxIndex; i++) {
      int sendId = macroElements.count(i);
540
541
542
543
      int recvId = 0;
      pdb.mpiComm.Allreduce(&sendId, &recvId, 1, MPI_INT, MPI_SUM);

      if (recvId != 1 && pdb.mpiRank == 0) {
544
545
546
	if (recvId == 0) {
	  ERROR_EXIT("Element %d has no member partition!\n", i);
	}
547

548
549
550
	if (recvId > 1) {
	  ERROR_EXIT("Element %d is member of more than pne partition!\n", i);
	}
551
552
553
554
555
      }
    }
  }


556
557
558
  void ParallelDebug::testDofContainerCommunication(MeshDistributor &pdb, 
						    RankToDofContainer &sendDofs,
						    RankToDofContainer &recvDofs)
559
  {
560
    FUNCNAME("ParallelDebug::testDofContainerCommunication()");
561
562
563
564
565
566
567
568
569

    std::map<int, int> sendNumber;
    for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it)
      sendNumber[it->first] = it->second.size();
    
    StdMpi<int> stdMpi(pdb.mpiComm);
    stdMpi.send(sendNumber);
    for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it)    
      stdMpi.recv(it->first);
570
    stdMpi.startCommunication();
571
572
     
    int foundError = 0;
573
    for (std::map<int, int>::iterator it = stdMpi.getRecvData().begin();
574
575
576
577
578
579
580
581
582
583
	 it != stdMpi.getRecvData().end(); ++it) {
      if (it->second != static_cast<int>(recvDofs[it->first].size())) {
	ERROR("Rank expectes %d DOFs to receive from rank %d, but %d DOFs are received!\n", 
	      recvDofs[it->first].size(), it->first, it->second);
	foundError = 1;
      }
    }

    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
584
585
586
  }


587
588
589
590
591
592
593
594
595
596
  void ParallelDebug::testDoubleDofs(Mesh *mesh)
  {
    FUNCNAME("ParallelDebug::testDoubleDofs()");

    std::map<WorldVector<double>, DegreeOfFreedom> cMap;
    int foundError = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
Thomas Witkowski's avatar
Merge    
Thomas Witkowski committed
597
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
598
599
600
601
602
	WorldVector<double> &c = elInfo->getCoord(i);
	if (cMap.count(c) == 0) {
	  cMap[c] = elInfo->getElement()->getDof(i, 0);
	} else {
	  if (cMap[c] != elInfo->getElement()->getDof(i, 0)) {
603
604
605
	    MSG("[DBG] Found two DOFs %d and %d with the same coords %f %f %f!\n",
		cMap[c], elInfo->getElement()->getDof(i, 0), 
		c[0], c[1], mesh->getDim() == 3 ? c[2] : 0.0);
606
607
608
609
610
611
612
613
614
615
616
617
618
	    foundError = 1;
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }
    
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
  }


619
  void ParallelDebug::printMapLocalGlobal(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
620
621
622
623
624
625
626
627
628
629
  {    
    if (rank == -1 || pdb.mpiRank == rank) {
      std::cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << std::endl;
      
      typedef std::map<DegreeOfFreedom, DegreeOfFreedom> DofMapping;
      typedef std::map<int, DofContainer> RankToDofContainer;

      for (DofMapping::iterator it = pdb.mapLocalGlobalDofs.begin();
	   it != pdb.mapLocalGlobalDofs.end(); it++) {
	DegreeOfFreedom localdof = -1;
630
631
	if (pdb.mapLocalDofIndex.count(it->first) > 0)
	  localdof = pdb.mapLocalDofIndex[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
	
	std::cout << "DOF " << it->first << " " 
		  << it->second << " " 
		  << localdof << std::endl;
	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(it->first, pdb.feSpace, coords);
	coords.print();

	for (RankToDofContainer::iterator rankit = pdb.sendDofs.begin();
	     rankit != pdb.sendDofs.end(); ++rankit) {
	  for (DofContainer::iterator dofit = rankit->second.begin();
	       dofit != rankit->second.end(); ++dofit)
	    if (**dofit == it->first)
	      std::cout << "SEND DOF TO " << rankit->first << std::endl;	  
	}

	for (RankToDofContainer::iterator rankit = pdb.recvDofs.begin();
	     rankit != pdb.recvDofs.end(); ++rankit) {
	  for (DofContainer::iterator dofit = rankit->second.begin();
	       dofit != rankit->second.end(); ++dofit)
	    if (**dofit == it->first)
	      std::cout << "RECV DOF FROM " << rankit->first << std::endl;	  
	}

	std::cout << "------" << std::endl;
      }
    }
  }


662
  void ParallelDebug::printMapPeriodic(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
663
  {
664
    FUNCNAME("ParallelDebug::printMapPeriodic()");
Thomas Witkowski's avatar
Thomas Witkowski committed
665

666
667
668
    ERROR_EXIT("Function must be rewritten!\n");

#if 0
Thomas Witkowski's avatar
Thomas Witkowski committed
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
    typedef std::map<DegreeOfFreedom, DegreeOfFreedom> DofMapping;
    typedef std::map<DegreeOfFreedom, std::set<DegreeOfFreedom> > PeriodicDofMap;

    if (rank == -1 || pdb.mpiRank == rank) {
      std::cout << "====== DOF MAP PERIODIC ====== " << std::endl;

      for (PeriodicDofMap::iterator it = pdb.periodicDof.begin();
	   it != pdb.periodicDof.end(); ++it) {
	std::cout << "DOF MAP " << it->first << ": ";
	for (std::set<DegreeOfFreedom>::iterator dofit = it->second.begin();
	     dofit != it->second.end(); ++dofit)
	  std::cout << *dofit << " ";
	std::cout << std::endl;

	DegreeOfFreedom localdof = -1;
	for (DofMapping::iterator dofIt = pdb.mapLocalGlobalDofs.begin();
	     dofIt != pdb.mapLocalGlobalDofs.end(); ++dofIt)
	  if (dofIt->second == it->first)
	    localdof = dofIt->first;

	TEST_EXIT(localdof != -1)("There is something wrong!\n");

	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(localdof, pdb.feSpace, coords);
	coords.print();
      }
    }
696
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
697
698
699
  }

  
700
701
702
703
  void ParallelDebug::printRankDofs(MeshDistributor &pdb, 
				    int rank, 
				    DofContainer& rankDofs,
				    DofContainer& rankAllDofs)
Thomas Witkowski's avatar
Thomas Witkowski committed
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
  {
    if (rank == -1 || pdb.mpiRank == rank) {
      std::cout << "====== RANK DOF INFORMATION ====== " << std::endl;

      std::cout << "  RANK OWNED DOFS: " << std::endl;
      for (DofContainer::iterator dofit = rankDofs.begin();
	   dofit != rankDofs.end(); ++dofit) {
	std::cout << "    " << **dofit << std::endl;
	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpace, coords);
	coords.print();
      }

      std::cout << "  RANK ALL DOFS: " << std::endl;
      for (DofContainer::iterator dofit = rankAllDofs.begin();
	   dofit != rankAllDofs.end(); ++dofit) {
	std::cout << "    " << **dofit << std::endl;
	WorldVector<double> coords;
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpace, coords);
	coords.print();
      }      
    }
  }

728

729
  void ParallelDebug::printBoundaryInfo(MeshDistributor &pdb)
730
  {
731
    FUNCNAME("ParallelDebug::printBoundaryInfo()");
732

Thomas Witkowski's avatar
Thomas Witkowski committed
733
    int tmp = 0;
734
    Parameters::get("parallel->debug->print boundary info", tmp);
Thomas Witkowski's avatar
Thomas Witkowski committed
735
    if (tmp <= 1)
736
737
      return;

738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
    for (InteriorBoundary::iterator it(pdb.myIntBoundary); !it.end(); ++it) {
      MSG("Rank owned boundary with rank %d: \n", it.getRank());
      MSG("  ranks obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->rankObj.elIndex, it->rankObj.subObj, it->rankObj.ithObj);
      MSG("  neigh obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
    }

    for (InteriorBoundary::iterator it(pdb.otherIntBoundary); !it.end(); ++it) {
      MSG("Other owned boundary with rank %d: \n", it.getRank());
      MSG("  ranks obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->rankObj.elIndex, it->rankObj.subObj, it->rankObj.ithObj);
      MSG("  neigh obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
    }
753
754

    for (InteriorBoundary::iterator it(pdb.periodicBoundary); !it.end(); ++it) {
755
756
      MSG("Periodic boundary (ID %d) with rank %d: \n", 
	  it->type, it.getRank());
757
758
759
760
761
      MSG("  ranks obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->rankObj.elIndex, it->rankObj.subObj, it->rankObj.ithObj);
      MSG("  neigh obj-ind: %d  sub-obj: %d   ith-obj: %d\n",
	  it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
    }    
762
763
  }

764

765
766
  void ParallelDebug::writeDebugFile(MeshDistributor &pdb,
				     std::string prefix, std::string postfix)
767
  {
768
    FUNCNAME("ParallelDebug::writeCoordsFile()");
769
770
771
772
773
774
775

    std::stringstream filename;
    filename << prefix << "-" << pdb.mpiRank << "." << postfix;

    DOFVector<WorldVector<double> > coords(pdb.feSpace, "tmp");
    pdb.mesh->getDofIndexCoords(pdb.feSpace, coords);

776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
    typedef std::map<int, std::vector<DegreeOfFreedom> > ElDofMap;
    ElDofMap elDofMap;
    TraverseStack stack;
    const BasisFunction *basisFcts = pdb.feSpace->getBasisFcts();
    std::vector<DegreeOfFreedom> localIndices(basisFcts->getNumber());
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      basisFcts->getLocalIndices(elInfo->getElement(), 
				 pdb.feSpace->getAdmin(), localIndices);
      elDofMap[elInfo->getElement()->getIndex()] = localIndices;
      elInfo = stack.traverseNext(elInfo);
    }

    // === Write informations about all DOFs. ===

791
792
    std::ofstream file;
    file.open(filename.str().c_str());
Thomas Witkowski's avatar
Thomas Witkowski committed
793
794
    file << "# First line contains number of DOFs, than each line has the format\n";
    file << "# Local DOF index     Global DOF index     Is rank DOF     x-coord     y-coord     z-coord\n";
795
796
797
    file << coords.getUsedSize() << "\n";
    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it) {
798
799
800
      file << it.getDOFIndex() << " " 
	   << pdb.mapLocalGlobalDofs[it.getDOFIndex()] << " "
	   << pdb.getIsRankDof(it.getDOFIndex());
801
802
803
804
      for (int i = 0; i < pdb.mesh->getDim(); i++)
	file << " " << (*it)[i];
      file << "\n";
    }
805
806
807

    // === Write to all elements in ranks mesh the included dofs. ===

808
809
810
    file << "\n\n";
    file << "# First line containes number of elements in mesh, second line contain the number of DOFs per element.\n";
    file << "# Than, each entry contains of two lines. The first is the element index, the second line is a list with the local DOF indices of this element.\n";
811
812
813
814
815
816
817
818
819
    file << elDofMap.size() << "\n";
    file << basisFcts->getNumber() << "\n";
    for (ElDofMap::iterator it = elDofMap.begin(); it != elDofMap.end(); ++it) {
      file << it->first << "\n";
      for (int i = 0; i < basisFcts->getNumber(); i++) 
	file << it->second[i] << " ";
      file << "\n";
    }

820
821
822
    file.close();
  }

823
824
825
826
827
828
829
830
831
832
833
834

  void ParallelDebug::writePartitioning(MeshDistributor &pdb, string filename)
  {
    FUNCNAME("ParallelDebug::writeParitioning()");

    map<int, double> vec;    
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    
    while (elInfo) {		  
      int index = elInfo->getElement()->getIndex();
835
      vec[index] = pdb.partitionMap[index];
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
      elInfo = stack.traverseNext(elInfo);
    }

    ElementFileWriter::writeFile(vec, pdb.mesh, filename);
  }

  void ParallelDebug::writePartitioningFile(std::string filename, 
					    int counter,
					    FiniteElemSpace *feSpace)
  {
    FUNCNAME("ParallelDebug::writePartitioningFile()");

    stringstream oss;
    oss << filename;
    if (counter >= 0)
      oss << "-" << counter;
    oss << ".vtu";

    DOFVector<double> tmpa(feSpace, "tmp");
    tmpa.set(MPI::COMM_WORLD.Get_rank());
    VtkWriter::writeFile(&tmpa, oss.str());
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
858
}