diff --git a/AMDiS/src/DualTraverse.cc b/AMDiS/src/DualTraverse.cc index 9a1e3bee26be40c7b7fdc81bdc4d687f41d70c90..28289844984123731b4fc7435f2d10a6490b6155 100644 --- a/AMDiS/src/DualTraverse.cc +++ b/AMDiS/src/DualTraverse.cc @@ -208,14 +208,9 @@ namespace AMDiS { if (largeElInfo->getLevel() == smallElInfo->getLevel()) return largeFace; - int refinementPathLength = smallElInfo->getRefinementPathLength(); - unsigned long refinementPath = smallElInfo->getRefinementPath(); - - TEST_EXIT_DBG(refinementPathLength != 0) + TEST_EXIT_DBG(smallElInfo->getRefinementPathLength() != 0) ("Refinement path is empty. This should not happen!\n"); - std::cout << "REFINEMENTPATHLENGTH = " << refinementPathLength << std::endl; - return -1; } } diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc index 250718e348ddacd3b36ae2eed835cec450ed982d..660741856fead04a5a439f7affc943e757b8f2a1 100644 --- a/AMDiS/src/Traverse.cc +++ b/AMDiS/src/Traverse.cc @@ -39,7 +39,7 @@ namespace AMDiS { traverse_fill_flag = fill_flag; TEST_EXIT_DBG(mesh)("No mesh!\n"); - TEST_EXIT_DBG(traverse_mesh->getMacroElements().size() > 0)("Mesh is empty!\n"); + TEST_EXIT(traverse_mesh->getMacroElements().size() > 0)("Mesh is empty!\n"); if (stack_size < 1) enlargeTraverseStack(); diff --git a/AMDiS/src/parallel/ElementObjectData.cc b/AMDiS/src/parallel/ElementObjectData.cc index 708574506697a05d82d2a9775b382f3bf308049e..075814aa0b7c055d8166596e9103257c9e8e25d3 100644 --- a/AMDiS/src/parallel/ElementObjectData.cc +++ b/AMDiS/src/parallel/ElementObjectData.cc @@ -422,10 +422,14 @@ namespace AMDiS { { FUNCNAME("ElementObjects::createReverseModeData()"); - // In 2D, all reverse modes are always true! + // === In 2D, all reverse modes are always true! === + if (mesh->getDim() == 2) return; + + // === First, create reverse modes for all "directly" neighbouring elements. === + for (map<DofEdge, vector<ElementObjectData> >::iterator edgeIt = edgeElements.begin(); edgeIt != edgeElements.end(); ++edgeIt) { @@ -473,6 +477,54 @@ namespace AMDiS { } } } + + + // === And create reverse modes for periodic neighbouring elements. === + + for (PerBoundMap<DofEdge>::iterator edgeIt = periodicEdges.begin(); + edgeIt != periodicEdges.end(); ++edgeIt) { + vector<ElementObjectData> &edges0 = edgeElements[edgeIt->first.first]; + vector<ElementObjectData> &edges1 = edgeElements[edgeIt->first.second]; + + for (unsigned int i = 0; i < edges0.size(); i++) { + BoundaryObject obj0(elIndexMap[edges0[i].elIndex], + elIndexTypeMap[edges0[i].elIndex], + EDGE, edges0[i].ithObject); + + for (unsigned int j = 0; j < edges1.size(); j++) { + BoundaryObject obj1(elIndexMap[edges1[j].elIndex], + elIndexTypeMap[edges1[j].elIndex], + EDGE, edges1[j].ithObject); + + bool reverseMode = + BoundaryObject::computeReverseMode(obj0, obj1, feSpace, edgeIt->second); + + edgeReverseMode[make_pair(edges0[i], edges1[j])] = reverseMode; + edgeReverseMode[make_pair(edges1[j], edges0[i])] = reverseMode; + } + } + } + + for (PerBoundMap<DofFace>::iterator faceIt = periodicFaces.begin(); + faceIt != periodicFaces.end(); ++faceIt) { + vector<ElementObjectData> &faces0 = faceElements[faceIt->first.first]; + vector<ElementObjectData> &faces1 = faceElements[faceIt->first.second]; + + TEST_EXIT_DBG(faces0.size() == faces1.size() == 1)("Should not happen!\n"); + + BoundaryObject obj0(elIndexMap[faces0[0].elIndex], + elIndexTypeMap[faces0[0].elIndex], + FACE, faces0[0].ithObject); + BoundaryObject obj1(elIndexMap[faces1[0].elIndex], + elIndexTypeMap[faces1[0].elIndex], + FACE, faces1[0].ithObject); + + bool reverseMode = + BoundaryObject::computeReverseMode(obj0, obj1, feSpace, faceIt->second); + + faceReverseMode[make_pair(faces0[0], faces1[0])] = reverseMode; + faceReverseMode[make_pair(faces1[0], faces0[0])] = reverseMode; + } } diff --git a/AMDiS/src/parallel/InteriorBoundary.cc b/AMDiS/src/parallel/InteriorBoundary.cc index 9b4366245b3c8f0dcc02b0b83966bdb1b01bfa9f..0e70b9e555b68474087ef3c2936aa3fafa5fbf72 100644 --- a/AMDiS/src/parallel/InteriorBoundary.cc +++ b/AMDiS/src/parallel/InteriorBoundary.cc @@ -92,16 +92,6 @@ namespace AMDiS { } - void BoundaryObject::setReverseMode(BoundaryObject &otherBound, - FiniteElemSpace *feSpace, - BoundaryType boundary) - { - FUNCNAME("BoundaryObject::setReverseMode()"); - - otherBound.reverseMode = computeReverseMode(*this, otherBound, feSpace, boundary); - } - - bool BoundaryObject::operator==(const BoundaryObject& other) const { return (other.elIndex == elIndex && diff --git a/AMDiS/src/parallel/InteriorBoundary.h b/AMDiS/src/parallel/InteriorBoundary.h index b68bff77c1722bb157827059fc7587a15aafbf78..a8b805b5aeada30a8e3ae56181786c5d5cab3b24 100644 --- a/AMDiS/src/parallel/InteriorBoundary.h +++ b/AMDiS/src/parallel/InteriorBoundary.h @@ -59,10 +59,6 @@ namespace AMDiS { FiniteElemSpace *feSpace, BoundaryType boundary); - void setReverseMode(BoundaryObject &otherBound, - FiniteElemSpace *feSpace, - BoundaryType boundary); - bool operator==(const BoundaryObject& other) const; bool operator!=(const BoundaryObject& other) const; diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index 2afa491d8ae2029358b707b767b14c0a2cd39f11..652ca9f0f0db52fdb7dae21474fc1727ef26c745 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -119,9 +119,9 @@ namespace AMDiS { TEST_EXIT(feSpace)("No FE space has been defined for the mesh distributor!\n"); TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n"); - int a; - char *b; - float zoltanVersion; + int a = 0; + char *b = NULL; + float zoltanVersion = 0.0; Zoltan_Initialize(a, &b, &zoltanVersion); elObjects.setFeSpace(feSpace); @@ -165,12 +165,9 @@ namespace AMDiS { setInitialElementWeights(); // and now partition the mesh - oldPartitionMap = partitioner->getPartitionMap(); - bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL); TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n"); - - partitionMap = partitioner->getPartitionMap(); + partitioner->createPartitionMap(partitionMap); #if (DEBUG != 0) @@ -742,6 +739,25 @@ namespace AMDiS { repartitionMesh(); nMeshChangesAfterLastRepartitioning = 0; } + + // === Print imbalance factor. === + + vector<int> nDofsInRank(mpiSize); + int nDofs = mesh->getDofAdmin(0).getUsedDofs(); + mpiComm.Gather(&nDofs, 1, MPI_INT, &(nDofsInRank[0]), 1, MPI_INT, 0); + + if (mpiRank == 0) { + int nOverallDofs = 0; + int maxDofs = numeric_limits<int>::min(); + for (int i = 0; i < mpiSize; i++) { + nOverallDofs += nDofsInRank[i]; + maxDofs = std::max(maxDofs, nDofsInRank[i]); + } + int avrgDofs = nOverallDofs / mpiSize; + double imbalance = (static_cast<double>(maxDofs - avrgDofs) / avrgDofs) * 100.0; + + MSG("Imbalancing factor: %.1f\%\n", imbalance); + } } @@ -935,7 +951,6 @@ namespace AMDiS { return; double timePoint = MPI::Wtime(); - repartitioningCounter++; #if (DEBUG != 0) ParallelDebug::testDoubleDofs(mesh); @@ -944,6 +959,7 @@ namespace AMDiS { ParallelDebug::writePartitioningFile(debugOutputDir + "partitioning", repartitioningCounter, feSpace); #endif + repartitioningCounter++; // === Create new element weights. === @@ -964,8 +980,6 @@ namespace AMDiS { MSG("Mesh partitioner created empty partition!\n"); return; } - oldPartitionMap = partitionMap; - // === Create map that maps macro element indices to pointers to the === // === macro elements. === @@ -1051,6 +1065,9 @@ namespace AMDiS { stdMpi.recv(it->first); stdMpi.startCommunication(); + if (interchangeVectors.size() == 0) + WARNING("There are no interchange vectors defined!\n"); + StdMpi<vector<vector<double> > > stdMpi2(mpiComm, true); stdMpi2.send(sendValues); for (std::map<int, std::vector<int> >::iterator it = partitioner->getRecvElements().begin(); @@ -1058,15 +1075,16 @@ namespace AMDiS { stdMpi2.recv(it->first); stdMpi2.startCommunication(); + // === Adapte the new macro elements due to the received mesh structure codes. === for (map<int, vector<int> >::iterator it = partitioner->getRecvElements().begin(); it != partitioner->getRecvElements().end(); ++it) { MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first]; vector<vector<double> > &recvValues = stdMpi2.getRecvData()[it->first]; - int i = 0; - int j = 0; - + int i = 0; + int j = 0; + TEST_EXIT_DBG(recvCodes.size() == it->second.size()) ("Should not happen!\n"); @@ -1132,17 +1150,15 @@ namespace AMDiS { // update the number of elements, vertices, etc. of the mesh. mesh->removeMacroElements(deleteMacroElements, feSpace); - // === Remove double DOFs. === MeshManipulation meshManipulation(feSpace); meshManipulation.deleteDoubleDofs(newMacroEl, elObjects); mesh->dofCompress(); - partitionMap = partitioner->getPartitionMap(); + partitioner->createPartitionMap(partitionMap); updateInteriorBoundaryInfo(); - updateLocalGlobalNumbering(); #if (DEBUG != 0) @@ -1215,11 +1231,11 @@ namespace AMDiS { } - elObjects.createReverseModeData(macroElIndexMap, macroElIndexTypeMap); - // === Create periodic data, if there are periodic boundary conditions. === elObjects.createPeriodicData(); + // === Create data about the reverse modes of neighbouring elements. === + elObjects.createReverseModeData(macroElIndexMap, macroElIndexTypeMap); // === Create mesh element data for this rank. === elObjects.createRankData(partitionMap); @@ -1283,7 +1299,6 @@ namespace AMDiS { AtomicBoundary& b = myIntBoundary.getNewAtomic(it2->first); b = bound; - // b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type); if (geoIndex == EDGE) b.neighObj.reverseMode = elObjects.getEdgeReverseMode(rankBoundEl, it2->second); if (geoIndex == FACE) @@ -1305,10 +1320,7 @@ namespace AMDiS { bound.type = INTERIOR; AtomicBoundary& b = otherIntBoundary.getNewAtomic(owner); - b = bound; - - // b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type); - + b = bound; if (geoIndex == EDGE) b.rankObj.reverseMode = elObjects.getEdgeReverseMode(rankBoundEl, ownerBoundEl); if (geoIndex == FACE) @@ -1391,9 +1403,9 @@ namespace AMDiS { b = bound; if (mpiRank > otherElementRank) - b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type); + b.neighObj.reverseMode = elObjects.getEdgeReverseMode(perEdgeEl0, perEdgeEl1); else - b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type); + b.rankObj.reverseMode = elObjects.getEdgeReverseMode(perEdgeEl0, perEdgeEl1); } } @@ -1408,26 +1420,26 @@ namespace AMDiS { TEST_EXIT_DBG(elObjects.getElements(it->first.second).size() == 1) ("Should not happen!\n"); - ElementObjectData& perEdgeEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank]; + ElementObjectData& perFaceEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank]; for (map<int, ElementObjectData>::iterator elIt = elObjects.getElementsInRank(it->first.second).begin(); elIt != elObjects.getElementsInRank(it->first.second).end(); ++elIt) { int otherElementRank = elIt->first; - ElementObjectData& perEdgeEl1 = elIt->second; + ElementObjectData& perFaceEl1 = elIt->second; AtomicBoundary bound; - bound.rankObj.el = macroElIndexMap[perEdgeEl0.elIndex]; - bound.rankObj.elIndex = perEdgeEl0.elIndex; - bound.rankObj.elType = macroElIndexTypeMap[perEdgeEl0.elIndex]; + bound.rankObj.el = macroElIndexMap[perFaceEl0.elIndex]; + bound.rankObj.elIndex = perFaceEl0.elIndex; + bound.rankObj.elType = macroElIndexTypeMap[perFaceEl0.elIndex]; bound.rankObj.subObj = FACE; - bound.rankObj.ithObj = perEdgeEl0.ithObject; + bound.rankObj.ithObj = perFaceEl0.ithObject; - bound.neighObj.el = macroElIndexMap[perEdgeEl1.elIndex]; - bound.neighObj.elIndex = perEdgeEl1.elIndex; - bound.neighObj.elType = macroElIndexTypeMap[perEdgeEl1.elIndex]; + bound.neighObj.el = macroElIndexMap[perFaceEl1.elIndex]; + bound.neighObj.elIndex = perFaceEl1.elIndex; + bound.neighObj.elType = macroElIndexTypeMap[perFaceEl1.elIndex]; bound.neighObj.subObj = FACE; - bound.neighObj.ithObj = perEdgeEl1.ithObject; + bound.neighObj.ithObj = perFaceEl1.ithObject; bound.type = it->second; @@ -1435,9 +1447,9 @@ namespace AMDiS { b = bound; if (mpiRank > otherElementRank) - b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type); + b.neighObj.reverseMode = elObjects.getFaceReverseMode(perFaceEl0, perFaceEl1); else - b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type); + b.rankObj.reverseMode = elObjects.getFaceReverseMode(perFaceEl0, perFaceEl1); } } @@ -1952,7 +1964,6 @@ namespace AMDiS { SerUtil::serialize(out, elemWeights); SerUtil::serialize(out, partitionMap); - SerUtil::serialize(out, oldPartitionMap); SerUtil::serialize(out, nRankDofs); SerUtil::serialize(out, nOverallDofs); @@ -1993,7 +2004,6 @@ namespace AMDiS { SerUtil::deserialize(in, elemWeights); SerUtil::deserialize(in, partitionMap); - SerUtil::deserialize(in, oldPartitionMap); SerUtil::deserialize(in, nRankDofs); SerUtil::deserialize(in, nOverallDofs); diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h index 1e760b569ca52e6179f3952e81db12fb8d1fce2a..2b14fced44413465c43146e1ababd21dcf9c98b4 100644 --- a/AMDiS/src/parallel/MeshDistributor.h +++ b/AMDiS/src/parallel/MeshDistributor.h @@ -31,6 +31,7 @@ #include "parallel/MeshPartitioner.h" #include "parallel/InteriorBoundary.h" +#include "AMDiS_fwd.h" #include "Global.h" #include "ProblemTimeInterface.h" #include "ProblemIterationInterface.h" @@ -38,7 +39,7 @@ #include "Serializer.h" #include "BoundaryManager.h" #include "ElementObjectData.h" -#include "AMDiS_fwd.h" +#include "SystemVector.h" namespace AMDiS { @@ -94,6 +95,13 @@ namespace AMDiS { interchangeVectors.push_back(vec); } + /// Adds all DOFVectors of a SystemVector to \ref interchangeVecs. + void addInterchangeVector(SystemVector *vec) + { + for (int i = 0; i < vec->getSize(); i++) + interchangeVectors.push_back(vec->getDOFVector(i)); + } + /** \brief * This function checks if the mesh has changed on at least on rank. In this case, * the interior boundaries are adapted on all ranks such that they fit together on @@ -448,12 +456,6 @@ namespace AMDiS { */ map<int, int> partitionMap; - /** \brief - * Stores an old partitioning of elements. To every macro element index the - * number of the rank it corresponds to is stored. - */ - map<int, int> oldPartitionMap; - /// Number of DOFs in the rank mesh. int nRankDofs; diff --git a/AMDiS/src/parallel/MeshPartitioner.h b/AMDiS/src/parallel/MeshPartitioner.h index 400df50d64ce5d8a0d1d1c3436d36b570b19aeeb..52701c7420c7f293520215cac03487d5519775f8 100644 --- a/AMDiS/src/parallel/MeshPartitioner.h +++ b/AMDiS/src/parallel/MeshPartitioner.h @@ -58,6 +58,8 @@ namespace AMDiS { virtual bool partition(map<int, double> &elemWeights, PartitionMode mode = INITIAL) = 0; + virtual void createPartitionMap(map<int, int>& partitionMap) = 0; + /// Write partitioner state to disk. void serialize(ostream &out); @@ -84,11 +86,6 @@ namespace AMDiS { return elementInRank; } - map<int, int>& getPartitionMap() - { - return partitionMap; - } - map<int, vector<int> >& getRecvElements() { return recvElements; diff --git a/AMDiS/src/parallel/ParMetisPartitioner.cc b/AMDiS/src/parallel/ParMetisPartitioner.cc index 4b65eeefe376d7b8c953522bd6aa4b9bef109e51..ee45e72c8f96d593ef98c71455d2f75a4d7e69e9 100644 --- a/AMDiS/src/parallel/ParMetisPartitioner.cc +++ b/AMDiS/src/parallel/ParMetisPartitioner.cc @@ -421,16 +421,11 @@ namespace AMDiS { // === Distribute new partition data. === - bool b = distributePartitioning(&(part[0])); - - if (b) - createPartitionMap(); - - return b; + return distributePartitioning(&(part[0])); } - void ParMetisPartitioner::createPartitionMap() + void ParMetisPartitioner::createPartitionMap(map<int, int>& pMap) { FUNCNAME("ParMetisPartitioner::createPartitionMap()"); @@ -468,6 +463,8 @@ namespace AMDiS { for (int i = 0; i < mpiSize; i++) for (int j = 0; j < nPartitionElements[i]; j++) partitionMap[partitionElements[elmdist[i] + j]] = i; + + pMap = partitionMap; } diff --git a/AMDiS/src/parallel/ParMetisPartitioner.h b/AMDiS/src/parallel/ParMetisPartitioner.h index 36db8e258c1e52e35c61393eb7a1276fda2f6825..cb1ef465ee366b5c9423d8167a4efdbbebcee23c 100644 --- a/AMDiS/src/parallel/ParMetisPartitioner.h +++ b/AMDiS/src/parallel/ParMetisPartitioner.h @@ -179,6 +179,8 @@ namespace AMDiS { bool partition(map<int, double> &elemWeights, PartitionMode mode = INITIAL); + void createPartitionMap(map<int, int>& partitionMap); + void setItr(float value) { itr = value; @@ -189,8 +191,6 @@ namespace AMDiS { // case, i.e., because there are empty partitions, the function returns false. bool distributePartitioning(int *part); - void createPartitionMap(); - protected: ParMetisMesh *parMetisMesh; diff --git a/AMDiS/src/parallel/ZoltanPartitioner.cc b/AMDiS/src/parallel/ZoltanPartitioner.cc index 7ebf62536d5243151a213f2c3576101ad0d8c0ea..c8c85fcb4ecbda12d0203ea9c338b4894236c262 100644 --- a/AMDiS/src/parallel/ZoltanPartitioner.cc +++ b/AMDiS/src/parallel/ZoltanPartitioner.cc @@ -54,9 +54,13 @@ namespace AMDiS { int *export_procs; int *export_to_part; - zoltan.Set_Param("LB_METHOD", "RCB"); - // zoltan.Set_Param("LB_METHOD", "GRAPH"); - // zoltan.Set_Param("LB_APPROACH", "PARTITION"); + if (mode == INITIAL) { + zoltan.Set_Param("LB_APPROACH", "PARTITION"); + } else { + zoltan.Set_Param("LB_APPROACH", "REPARTITION"); + zoltan.Set_Param("LB_METHOD", "GRAPH"); + } + zoltan.Set_Param("OBJ_WEIGHT_DIM", "1"); #if (DEBUG != 0) diff --git a/AMDiS/src/parallel/ZoltanPartitioner.h b/AMDiS/src/parallel/ZoltanPartitioner.h index 41a6319fd598e0413fc8ef2371f5f79ca5725284..909c2c8fdacab4e49eb47f4d83c1f4754111fab0 100644 --- a/AMDiS/src/parallel/ZoltanPartitioner.h +++ b/AMDiS/src/parallel/ZoltanPartitioner.h @@ -42,6 +42,11 @@ namespace AMDiS { /// \ref MeshPartitioner::partition bool partition(map<int, double> &weights, PartitionMode mode = INITIAL); + void createPartitionMap(map<int, int> &pMap) + { + pMap = partitionMap; + } + protected: void createPartitionMap();