ParallelDebug.cc 32 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"
16
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
17
18
19
20
#include "DOFVector.h"
#include "FixVec.h"
#include "StdMpi.h"
#include "Debug.h"
21
#include "io/VtkWriter.h"
22
#include "ElementDofIterator.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
23
24
25

namespace AMDiS {

26
27
28
  using namespace std;


29
  void ParallelDebug::testInteriorBoundary(MeshDistributor &pdb)
Thomas Witkowski's avatar
Thomas Witkowski committed
30
  {
31
    FUNCNAME("ParallelDebug::testInteriorBoundary()");
Thomas Witkowski's avatar
Thomas Witkowski committed
32

33
    vector<int*> sendBuffers, recvBuffers;
Thomas Witkowski's avatar
Thomas Witkowski committed
34

35
36
37
    MPI::Request request[pdb.intBoundary.own.size() + 
			 pdb.intBoundary.other.size() +
                         pdb.intBoundary.periodic.size() * 2];
Thomas Witkowski's avatar
Thomas Witkowski committed
38
39
    int requestCounter = 0;

40
41
42

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

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

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

57
58
59

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

60
61
    for (RankToBoundMap::iterator rankIt = pdb.intBoundary.other.begin();
	 rankIt != pdb.intBoundary.other.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
62
63
64
65
66
67
68
69
70
      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);
    }


71
72
    // === To the last, do the same of periodic boundaries. ===

73
74
    for (RankToBoundMap::iterator rankIt = pdb.intBoundary.periodic.begin();       
	 rankIt != pdb.intBoundary.periodic.end(); ++rankIt) {
75
76
77
      if (rankIt->first == pdb.mpiRank)
	continue;

78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
      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
98
99
100
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i];

101
102
103
104
105

    // === 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
106
    int bufCounter = 0;
107
108
    for (RankToBoundMap::iterator rankIt = pdb.intBoundary.other.begin();
	 rankIt != pdb.intBoundary.other.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
109

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

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

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

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


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

130
131
      for (unsigned int i = 0; i < rankIt->second.size(); i++) {
	int elIndex1 = recvBuffers[bufCounter][i];
132
	int elIndex2 = pdb.intBoundary.periodic[rankIt->first][i].neighObj.elIndex;
133

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

      delete [] recvBuffers[bufCounter++];
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
142
143
144
  }


145
146
147
148
  void ParallelDebug::testPeriodicBoundary(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testPeriodicBoundary()");

Thomas Witkowski's avatar
Thomas Witkowski committed
149
150
151
    MSG(" NO TEST HERE!\n");
//     for (unsigned int i = 0; i < pdb.feSpaces.size(); i++)
//       testPeriodicBoundary(pdb, pdb.feSpaces[i]);
152
153
  }

154

155
156
157
158
  void ParallelDebug::testPeriodicBoundary(MeshDistributor &pdb,
					   const FiniteElemSpace *feSpace)
  {
    FUNCNAME("ParallelDebug::testPeriodicBoundary()");
159

160
161
    // === 1. check: All periodic DOFs should have at least a correct number ===
    // === of periodic associations.                                         ===
162
163
164
165
166

    PeriodicMap &perMap = pdb.getPeriodicMap();
    for (map<int, std::set<BoundaryType> >::iterator it = 
	   perMap.periodicDofAssociations[feSpace].begin();
	 it != perMap.periodicDofAssociations[feSpace].end(); ++it) {
167
      WorldVector<double> c;
168
      pdb.mesh->getDofIndexCoords(it->first, pdb.feSpaces[0], c);
169
170
171
172
173
174
175
      int nAssoc = it->second.size();
    }    


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

176
177
178
179
180
181
    StdMpi<PeriodicDofMap> stdMpi(pdb.mpiComm, true);

    if (pdb.mpiRank == 0) {
      for (int i = 1; i < pdb.mpiSize; i++)
	stdMpi.recv(i);
    } else {
182
      stdMpi.send(0, perMap.periodicDofMap[feSpace]);
183
184
    }

185
    stdMpi.startCommunication();
186
187
188
189
190
191

    int foundError = 0;

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

    if (pdb.mpiRank == 0) {
192
      // Stores to each rank the periodic DOF mappings of this rank.
193
      map<int, PeriodicDofMap> rankToMaps;
194
      PeriodicDofMap dofMap = perMap.periodicDofMap[feSpace];
195
196
197
198
199
200
201
202
      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) {
Thomas Witkowski's avatar
Thomas Witkowski committed
203
	  for (std::map<DegreeOfFreedom, DegreeOfFreedom>::iterator dofIt = it->second.begin();
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
	       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) {
Thomas Witkowski's avatar
Thomas Witkowski committed
222
	for (std::map<DegreeOfFreedom, DegreeOfFreedom>::iterator dofIt = it->second.begin();
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
	     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");
251
252


253
254
255
256
    // === 3. check: On all edge and face periodic DOFs, at least one        ===
    // === componant of coordinates of each periodic DOF pair must be equal  ===
    // === (at least as long we consider periodic boundaries only on         ===
    // === rectangulars and boxes.                                           ===
257
258

    RankToCoords sendCoords;
259
    map<int, vector<BoundaryType> > rankToDofType;
260

261
262
    for (RankToBoundMap::iterator it = pdb.intBoundary.periodic.begin();
	 it != pdb.intBoundary.periodic.end(); ++it) {
263
264
265
266
267
268
269
270
271
      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;
272
	boundIt->rankObj.el->getAllDofs(feSpace, boundIt->rankObj, dofs);
273
274
275
	
	for (unsigned int i = 0; i < dofs.size(); i++) {
	  WorldVector<double> c;
276
	  pdb.mesh->getDofIndexCoords(*(dofs[i]), feSpace, c);
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
	  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);   
294
    stdMpiCoords.startCommunication();
295
296
297
298
299
300
301
302
303
304
305
306
307
308


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

309
	if (nEqual == 0) {
310
311
312
313
314
315
316
317
318
319
320
321
322
323
	  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");
324
325
326
  }


327
  void ParallelDebug::testCommonDofs(MeshDistributor &pdb, bool printCoords)
Thomas Witkowski's avatar
Thomas Witkowski committed
328
  {
329
    FUNCNAME("ParallelDebug::testCommonDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331

    clock_t first = clock();
332
    // Get FE space with basis functions of the highest degree
333
    const FiniteElemSpace *feSpace = pdb.feSpaces[pdb.feSpaces.size() - 1];
Thomas Witkowski's avatar
Thomas Witkowski committed
334
335

    int testCommonDofs = 1;
336
    Parameters::get("dbg->test common dofs", testCommonDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
337
338
339
340
341
342
    if (testCommonDofs == 0) {
      MSG("Skip test common dofs!\n");
      return;
    }

    /// Defines a mapping type from rank numbers to sets of DOFs.
343
    typedef map<int, DofContainer> RankToDofContainer;
Thomas Witkowski's avatar
Thomas Witkowski committed
344
345

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

351
352
    // 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
353
354
355
    // this rank expectes to receive from.
    RankToCoords recvCoords;

356
    DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds");
Thomas Witkowski's avatar
Thomas Witkowski committed
357
    pdb.mesh->getDofIndexCoords(coords);
358

359
360
    for (DofComm::Iterator it(pdb.dofComm.getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
361
362
      for (; !it.endDofIter(); it.nextDof())
	sendCoords[it.getRank()].push_back(coords[it.getDofIndex()]);
Thomas Witkowski's avatar
Thomas Witkowski committed
363

364
365
    for (DofComm::Iterator it(pdb.dofComm.getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank())
366
367
      for (; !it.endDofIter(); it.nextDof())
	recvCoords[it.getRank()].push_back(coords[it.getDofIndex()]);
Thomas Witkowski's avatar
Thomas Witkowski committed
368

369
370
371
    vector<int> sendSize(pdb.mpiSize, 0);
    vector<int> recvSize(pdb.mpiSize, 0);
    vector<int> recvSizeBuffer(pdb.mpiSize, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
374
375
376
377
378
379
380
381
382
383
384
    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;

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

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

393
394
      request[requestCounter++] = 
	pdb.mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
395
396
397
398
399
    }

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


400
    int foundError = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
401
402
403
404
405
    for (int i = 0; i < pdb.mpiSize; i++) {
      if (i == pdb.mpiRank)
	continue;

      if (recvSize[i] != recvSizeBuffer[i]) {
406
407
408
	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
409
410
      }
    }
411
412
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
415
416
417
418
419
420

    // === 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);   
421
    stdMpi.startCommunication();
Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
424
425
426
427
428

    int dimOfWorld = Global::getGeo(WORLD);

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

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

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

466
    MSG("Test common dofs needed %.5f seconds\n", TIME_USED(first, clock()));
Thomas Witkowski's avatar
Thomas Witkowski committed
467
468
469
  }


470
471
472
473
  void ParallelDebug::testGlobalIndexByCoords(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::testGlobalIndexByCoords()");

474
    // Get FE space with basis functions of the highest degree
475
    const FiniteElemSpace *feSpace = pdb.feSpaces[pdb.feSpaces.size() - 1];
476
477

    DOFVector<WorldVector<double> > coords(feSpace, "tmp");
Thomas Witkowski's avatar
Thomas Witkowski committed
478
    pdb.mesh->getDofIndexCoords(coords);
479

480
    typedef map<WorldVector<double>, int> CoordsIndexMap;
481
482
483
    CoordsIndexMap coordsToIndex;

    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
484
    for (it.reset(); !it.end(); ++it) {
485
      coordsToIndex[(*it)] = (*(pdb.dofMaps[0]))[feSpace][it.getDOFIndex()].global;
486
487
488
489
      // MSG("   CHECK FOR DOF %d/%d AT COORDS %f %f %f\n",
      // 	  it.getDOFIndex(),
      // 	  coordsToIndex[(*it)], 
      // 	  (*it)[0], (*it)[1], (pdb.mesh->getDim() == 3 ? (*it)[2] : 0.0));
490
491
    }

492
    StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true);
493
494
    for (DofComm::Iterator it(pdb.dofComm.getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
495
      stdMpi.send(it.getRank(), coordsToIndex);
496
497
    for (DofComm::Iterator it(pdb.dofComm.getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank())
498
      stdMpi.recv(it.getRank());
499
   
500
    stdMpi.startCommunication();
501
502

    int foundError = 0;
503
504
    for (DofComm::Iterator it(pdb.dofComm.getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
505
      CoordsIndexMap& otherCoords = stdMpi.getRecvData(it.getRank());
506
507
508
509
510

      for (CoordsIndexMap::iterator coordsIt = otherCoords.begin();
	   coordsIt != otherCoords.end(); ++coordsIt) {
	if (coordsToIndex.count(coordsIt->first) == 1 &&
	    coordsToIndex[coordsIt->first] != coordsIt->second) {
511
	  stringstream oss;
512
513
514
515
516
517
	  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: " 
518
	      << coordsToIndex[coordsIt->first] << ") and on rank "
519
	      << it.getRank() << " (global index: " << coordsIt->second << ")";
520
521
522
523
524
525
526
527
528
529
530

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

531

532
  void ParallelDebug::testAllElements(MeshDistributor &pdb)
533
  {
534
    FUNCNAME("ParallelDebug::testAllElements()");
535

536
    std::set<int> macroElements;
537
538
    int minElementIndex = numeric_limits<int>::max();
    int maxElementIndex = numeric_limits<int>::min();   
539

540
    TraverseStack stack;
541
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, 0, Mesh::CALL_EL_LEVEL);
542
    while (elInfo) {
543
544
545
546
547
548
549
550
551
552
      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);
553

554
555
556
    TEST_EXIT(globalMinIndex == 0)("No macro element with index 0!\n");
    for (int i = 0; i <= globalMaxIndex; i++) {
      int sendId = macroElements.count(i);
557
558
559
560
      int recvId = 0;
      pdb.mpiComm.Allreduce(&sendId, &recvId, 1, MPI_INT, MPI_SUM);

      if (recvId != 1 && pdb.mpiRank == 0) {
561
562
563
	if (recvId == 0) {
	  ERROR_EXIT("Element %d has no member partition!\n", i);
	}
564

565
566
567
	if (recvId > 1) {
	  ERROR_EXIT("Element %d is member of more than pne partition!\n", i);
	}
568
569
570
571
572
      }
    }
  }


573
  void ParallelDebug::testDofContainerCommunication(MeshDistributor &pdb)
574
  {
575
    FUNCNAME("ParallelDebug::testDofContainerCommunication()");    
576

577
578
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;

579
    map<int, int> sendNumber;
580
581
582
583
    for (it_type it = pdb.dofComm.getSendDofs()[0].begin(); 
	 it != pdb.dofComm.getSendDofs()[0].end(); ++it)
      for (map<const FiniteElemSpace*, DofContainer>::iterator dcIt = it->second.begin(); 
	   dcIt != it->second.end(); ++dcIt)
584
	sendNumber[it->first] += dcIt->second.size();
585
       
586
    map<int, int> recvNumber;
587
588
589
590
    for (it_type it = pdb.dofComm.getRecvDofs()[0].begin(); 
	 it != pdb.dofComm.getRecvDofs()[0].end(); ++it)
      for (map<const FiniteElemSpace*, DofContainer>::iterator dcIt = it->second.begin(); 
	   dcIt != it->second.end(); ++dcIt)
591
	recvNumber[it->first] += dcIt->second.size();
592
593
594
    
    StdMpi<int> stdMpi(pdb.mpiComm);
    stdMpi.send(sendNumber);
595
596
    for (it_type it = pdb.dofComm.getRecvDofs()[0].begin(); 
	 it != pdb.dofComm.getRecvDofs()[0].end(); ++it)
597
      stdMpi.recv(it->first);
598
    stdMpi.startCommunication();
599

600
    int foundError = 0;
601
    for (map<int, int>::iterator it = stdMpi.getRecvData().begin();
602
	 it != stdMpi.getRecvData().end(); ++it) {
603
      if (it->second != recvNumber[it->first]) {
604
	ERROR("Rank expectes %d DOFs to receive from rank %d, but %d DOFs are received!\n", 
605
	      recvNumber[it->first], it->first, it->second);
606
607
608
609
610
611
	foundError = 1;
      }
    }

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


615
616
617
618
  void ParallelDebug::testDoubleDofs(Mesh *mesh)
  {
    FUNCNAME("ParallelDebug::testDoubleDofs()");

619
    map<WorldVector<double>, DegreeOfFreedom> cMap;
620
621
622
623
624
    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
625
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
626
627
628
629
630
	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)) {
631
632
633
	    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);
634
635
636
637
638
639
640
641
642
643
644
645
646
	    foundError = 1;
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }
    
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
  }


647
  void ParallelDebug::printMapLocalGlobal(MeshDistributor &pdb, int rank)
648
649
650
651
652
653
  {   
    FUNCNAME("ParallelDebug::printMapLocalGlobal()");

    ERROR_EXIT("Rewrite this function!\n");

#if 0
Thomas Witkowski's avatar
Thomas Witkowski committed
654
    if (rank == -1 || pdb.mpiRank == rank) {
655
      const FiniteElemSpace *feSpace = pdb.feSpaces[0];
656

657
      cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << endl;
658
            
659
660
      DofMap &dofMap = pdb.dofMap[feSpace].getMap();
      for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
661
	cout << "DOF " << it->first << " " << it->second.global << "\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
662
	WorldVector<double> coords;
663
	pdb.mesh->getDofIndexCoords(it->first, feSpace, coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
	coords.print();

666
	for (DofComm::Iterator rit(pdb.dofComm.getSendDofs(), feSpace); 
667
	     !rit.end(); rit.nextRank())
668
669
670
671
	  for (; !rit.endDofIter(); rit.nextDof())
	    if (it->first == rit.getDofIndex())
	      cout << "SEND DOF TO " << rit.getRank() << endl;
	
672
	for (DofComm::Iterator rit(pdb.dofComm.getRecvDofs(), feSpace); 
673
	     !rit.end(); rit.nextRank())
674
675
676
	  for (; !rit.endDofIter(); rit.nextDof())
	    if (it->first == rit.getDofIndex())
	      cout << "RECV DOF FROM " << rit.getRank() << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
677

678
	cout << "------" << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680
      }
    }
681
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
682
683
684
  }


685
  void ParallelDebug::printMapPeriodic(MeshDistributor &pdb, int rank)
Thomas Witkowski's avatar
Thomas Witkowski committed
686
  {
687
    FUNCNAME("ParallelDebug::printMapPeriodic()");
Thomas Witkowski's avatar
Thomas Witkowski committed
688

Thomas Witkowski's avatar
Thomas Witkowski committed
689
    ERROR_EXIT("Function must be rewritten, check svn for old code!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
690
691
692
  }

  
693
694
695
696
  void ParallelDebug::printRankDofs(MeshDistributor &pdb, 
				    int rank, 
				    DofContainer& rankDofs,
				    DofContainer& rankAllDofs)
Thomas Witkowski's avatar
Thomas Witkowski committed
697
698
  {
    if (rank == -1 || pdb.mpiRank == rank) {
699
      cout << "====== RANK DOF INFORMATION ====== " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
700

701
      cout << "  RANK OWNED DOFS: " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
702
703
      for (DofContainer::iterator dofit = rankDofs.begin();
	   dofit != rankDofs.end(); ++dofit) {
704
	cout << "    " << **dofit << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
705
	WorldVector<double> coords;
706
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpaces[0], coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
707
708
709
	coords.print();
      }

710
      cout << "  RANK ALL DOFS: " << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
711
712
      for (DofContainer::iterator dofit = rankAllDofs.begin();
	   dofit != rankAllDofs.end(); ++dofit) {
713
	cout << "    " << **dofit << endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
714
	WorldVector<double> coords;
715
	pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpaces[0], coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
718
719
720
	coords.print();
      }      
    }
  }

721

Thomas Witkowski's avatar
Thomas Witkowski committed
722
  void ParallelDebug::printBoundaryInfo(InteriorBoundary &intBoundary,
723
724
					int level, 
					bool force)
725
  {
726
    FUNCNAME("ParallelDebug::printBoundaryInfo()");
727

Thomas Witkowski's avatar
Thomas Witkowski committed
728
    int tmp = 0;
729
    Parameters::get("parallel->debug->print boundary info", tmp);
730
    if (tmp <= 0 && force == false)
731
732
      return;

733
734
    MSG("Interior boundary info:\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
735
    for (InteriorBoundary::iterator it(intBoundary.own, level); 
736
	 !it.end(); ++it) {
737
738
739
740
741
742
743
      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);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
744
    for (InteriorBoundary::iterator it(intBoundary.other, level); 
745
	 !it.end(); ++it) {
746
747
748
749
750
751
      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);
    }
752

Thomas Witkowski's avatar
Thomas Witkowski committed
753
    for (InteriorBoundary::iterator it(intBoundary.periodic, level); 
754
	 !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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
765
766
767
768
  void ParallelDebug::writeDebugFile(const FiniteElemSpace *feSpace,
				     ParallelDofMapping &dofMap,
				     string prefix, 
				     string postfix)
769
  {
770
    FUNCNAME("ParallelDebug::writeCoordsFile()");
771

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
772
    Mesh *mesh = feSpace->getMesh();
773

774
    stringstream filename;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
775
    filename << prefix << "-" << MPI::COMM_WORLD.Get_rank() << "." << postfix;
776

777
    DOFVector<WorldVector<double> > coords(feSpace, "tmp");
Thomas Witkowski's avatar
Thomas Witkowski committed
778
    mesh->getDofIndexCoords(coords);
779

780
    typedef map<int, vector<DegreeOfFreedom> > ElDofMap;
781
782
    ElDofMap elDofMap;
    TraverseStack stack;
783
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
784
    vector<DegreeOfFreedom> localIndices(basisFcts->getNumber());
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
785
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
786
787
    while (elInfo) {
      basisFcts->getLocalIndices(elInfo->getElement(), 
788
				 feSpace->getAdmin(), localIndices);
789
790
791
792
793
794
      elDofMap[elInfo->getElement()->getIndex()] = localIndices;
      elInfo = stack.traverseNext(elInfo);
    }

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

795
    ofstream file;
796
    file.open(filename.str().c_str());
Thomas Witkowski's avatar
Thomas Witkowski committed
797
//     file << "# First line contains number of DOFs, than each line has the format\n";
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
798
799
//     file << "# DOF index     Local DOF index     Global DOF index     Is rank DOF     x-coord     y-coord     z-coord\n";
    file << dofMap[feSpace].size() << "\n";
800
801
    DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
    for (it.reset(); !it.end(); ++it) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
802
803
804
805
806
807
808
809
810
      if (dofMap[feSpace].isSet(it.getDOFIndex())) {
	file << it.getDOFIndex() << " " 
	     << dofMap[feSpace][it.getDOFIndex()].local << " "
	     << dofMap[feSpace][it.getDOFIndex()].global << " "
	     << dofMap[feSpace].isRankDof(it.getDOFIndex());
	for (int i = 0; i < mesh->getDim(); i++)
	  file << " " << (*it)[i];
	file << "\n";
      }
811
    }
812
813
814

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

Thomas Witkowski's avatar
Thomas Witkowski committed
815
816
817
//     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";
818
819
820
821
822
823
824
825
826
    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";
    }

827
828
829
    file.close();
  }

830
831
832
833
834
835
836

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

    map<int, double> vec;    
    TraverseStack stack;
837
    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, Mesh::CALL_LEAF_EL);
838
839
840
    
    while (elInfo) {		  
      int index = elInfo->getElement()->getIndex();
841
      vec[index] = pdb.partitionMap[index];
842
843
844
845
846
847
      elInfo = stack.traverseNext(elInfo);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
848

849
  void ParallelDebug::writePartitioningFile(string filename, 
850
					    int counter,
851
					    const FiniteElemSpace *feSpace)
852
853
854
855
856
857
858
859
860
  {
    FUNCNAME("ParallelDebug::writePartitioningFile()");

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

861
862
863
    DOFVector<double> tmp(feSpace, "tmp");
    tmp.set(MPI::COMM_WORLD.Get_rank());
    VtkWriter::writeFile(&tmp, oss.str());
864
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886

  
  bool ParallelDebug::followThisBound(int rankElIndex, int neighElIndex)
  {
    FUNCNAME("ParallelDebug::followThisBound()");

    int el0 = std::min(rankElIndex, neighElIndex);
    int el1 = std::max(rankElIndex, neighElIndex);

    vector<int> els;
    Parameters::get("parallel->debug->follow boundary", els);
    if (els.size() != 2)
      return false;

    return (el0 == els[0] && el1 == els[1]);
  }


  void ParallelDebug::followBoundary(MeshDistributor &pdb)
  {
    FUNCNAME("ParallelDebug::followBoundary()");

887
    for (InteriorBoundary::iterator it(pdb.intBoundary.own); !it.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
888
889
890
      if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex))
	debug::writeLocalElementDofs(pdb.mpiRank, 
				     it->rankObj.elIndex, 
891
				     pdb.feSpaces[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
892

893
    for (InteriorBoundary::iterator it(pdb.intBoundary.other); !it.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
894
895
896
      if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex))
	debug::writeLocalElementDofs(pdb.mpiRank,
				     it->rankObj.elIndex,
897
				     pdb.feSpaces[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
  }


  void ParallelDebug::followBoundary(Mesh *mesh, 
				     AtomicBoundary &bound, 
				     MeshStructure &code)
  {
    FUNCNAME("ParallelDebug::followBoundary()");

    if (mesh->getDim() != bound.rankObj.subObj)
      return;

    if (!followThisBound(bound.rankObj.elIndex, bound.neighObj.elIndex))
      return;

    MSG("Mesh structure code of bound %d/%d <-> %d/%d: %s\n",
	bound.rankObj.elIndex, mesh->getDim(),
	bound.neighObj.elIndex, mesh->getDim(), 
	code.toStr().c_str());
  }
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
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

  void ParallelDebug::writeCsvElementMap(const FiniteElemSpace *feSpace,
			       ParallelDofMapping &dofMap,
			       string prefix, 
			       string postfix)
  {
	FUNCNAME("ParallelDebug::writeCsvElementMap()");

	MSG("writing local Element map to CSV File \n");

	Mesh *mesh = feSpace->getMesh();

	stringstream filename;
    filename << prefix << "-" << MPI::COMM_WORLD.Get_rank() << "." << postfix;

	ofstream file;
    file.open(filename.str().c_str());

	DOFVector<WorldVector<double> > dofCoords(feSpace, "DOF coords");
    mesh->getDofIndexCoords(dofCoords);

	TraverseStack stack;
	ElInfo *elInfo = stack.traverseFirst(mesh,	
									-1, 
									Mesh::CALL_EVERY_EL_POSTORDER);

	while (elInfo) {
		/*
		MSG("Start traverse element %d in level %d\n",
			elInfo->getElement()->getIndex(),
			elInfo->getLevel());
		 */

		// check if Element is Leaflevel
		// because higherOrderDofs(=NON VERTICES DOFS) are not saved in nonleaf Elements
		if (elInfo->getElement()->isLeaf()) {

			//traverse ALL DOFS
  			ElementDofIterator elDofIter(feSpace, true);
			elDofIter.reset(elInfo->getElement());
	
			int locDofNumber = 0;  

			do {
   	
		    	WorldVector<double> c = dofCoords[elDofIter.getDof()];
	
				file << elInfo->getElement()->getIndex() << "\t";	//element number
				file << elInfo->getLevel() << "\t";					//element Level
				file << elDofIter.getDof() << "\t";					//dof number
				file << elDofIter.getCurrentPos()+1 << "\t";		//dof type (+1 to pseudo convert to geoIndex, does not work for CENTER DOFS( geoIndex=0))
				file << c[0] << "\t";								//dof coords
				file << c[1] << "\t";	
				file << c[2] << "\t";
				file << locDofNumber << "\t";						//local Dof number
				file << elInfo->getType() << "\n";					//element Type
		
				locDofNumber++;
		    } while (elDofIter.next());
		} else {
	
			// traverse only VERTEX DOFS
			for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
				DegreeOfFreedom dof = elInfo->getElement()->getDof(i, 0);	

				WorldVector<double> c = dofCoords[dof];
		
				file << elInfo->getElement()->getIndex() << "\t";	//element number
				file << elInfo->getLevel() << "\t";					//element Level
				file << dof << "\t";								//dof number
				file << VERTEX << "\t";								//dof type
				file << c[0] << "\t";								//dof coords
				file << c[1] << "\t";	
				file << c[2] << "\t";
				file << i << "\t";									//local Dof number
				file << elInfo->getType() << "\n";					//element Type
			}		
		}

   		elInfo = stack.traverseNext(elInfo);
	}

  }
For faster browsing, not all history is shown. View entire blame