1 #ifndef DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH 2 #define DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH 3 4 /** \file 5 * \brief Compute a repartitioning of a Dune grid using ParMetis 6 */ 7 8 #include <algorithm> 9 #include <vector> 10 11 #include <dune/common/parallel/mpihelper.hh> 12 #include <dune/common/exceptions.hh> 13 14 #include <dune/geometry/referenceelements.hh> 15 16 #include <dune/grid/utility/globalindexset.hh> 17 #include <dune/grid/common/mcmgmapper.hh> 18 #include <dune/grid/common/rangegenerators.hh> 19 20 #if HAVE_PARMETIS 21 22 #include <parmetis.h> 23 24 // only enable for ParMETIS because the implementation uses functions that 25 // are not emulated by scotch 26 #ifdef PARMETIS_MAJOR_VERSION 27 28 namespace Dune 29 { 30 31 /** \brief Compute a repartitioning of a Dune grid using ParMetis 32 * 33 * \author Benjamin Bykowski 34 */ 35 template<class GridView> 36 struct ParMetisGridPartitioner { 37 38 // define index type as provided by ParMETIS 39 #if PARMETIS_MAJOR_VERSION > 3 40 typedef idx_t idx_type; 41 typedef ::real_t real_type; 42 #else 43 typedef int idx_type; 44 typedef float real_type; 45 #endif // PARMETIS_MAJOR_VERSION > 3 46 47 enum { 48 dimension = GridView::dimension 49 }; 50 51 52 /** \brief Create an initial partitioning of a Dune grid, i.e., not taking into account communication cost 53 * 54 * This pipes a Dune grid into the method ParMETIS_V3_PartMeshKway (see the ParMetis documentation for details) 55 * 56 * \param gv The grid view to be partitioned 57 * \param mpihelper The MPIHelper object, needed to get the MPI communicator 58 * 59 * \return std::vector with one uint per All_Partition element. For each element, the entry is the 60 * number of the partition the element is assigned to. 61 */ partitionDune::ParMetisGridPartitioner62 static std::vector<unsigned> partition(const GridView& gv, const Dune::MPIHelper& mpihelper) { 63 const unsigned numElements = gv.size(0); 64 65 std::vector<unsigned> part(numElements); 66 67 // Setup parameters for ParMETIS 68 idx_type wgtflag = 0; // we don't use weights 69 idx_type numflag = 0; // we are using C-style arrays 70 idx_type ncon = 1; // number of balance constraints 71 idx_type ncommonnodes = 2; // number of nodes elements must have in common to be considered adjacent to each other 72 idx_type options[4] = {0, 0, 0, 0}; // use default values for random seed, output and coupling 73 idx_type edgecut; // will store number of edges cut by partition 74 idx_type nparts = mpihelper.size(); // number of parts equals number of processes 75 std::vector<real_type> tpwgts(ncon*nparts, 1./nparts); // load per subdomain and weight (same load on every process) 76 std::vector<real_type> ubvec(ncon, 1.05); // weight tolerance (same weight tolerance for every weight there is) 77 78 // The difference elmdist[i+1] - elmdist[i] is the number of nodes that are on process i 79 std::vector<idx_type> elmdist(nparts+1); 80 elmdist[0] = 0; 81 std::fill(elmdist.begin()+1, elmdist.end(), gv.size(0)); // all elements are on process zero 82 83 // Create and fill arrays "eptr", where eptr[i] is the number of vertices that belong to the i-th element, and 84 // "eind" contains the vertex-numbers of the i-the element in eind[eptr[i]] to eind[eptr[i+1]-1] 85 std::vector<idx_type> eptr, eind; 86 int numVertices = 0; 87 eptr.push_back(numVertices); 88 89 for (const auto& element : elements(gv, Partitions::interior)) { 90 const size_t curNumVertices = referenceElement<double, dimension>(element.type()).size(dimension); 91 92 numVertices += curNumVertices; 93 eptr.push_back(numVertices); 94 95 for (size_t k = 0; k < curNumVertices; ++k) 96 eind.push_back(gv.indexSet().subIndex(element, k, dimension)); 97 } 98 99 // Partition mesh using ParMETIS 100 if (0 == mpihelper.rank()) { 101 MPI_Comm comm = Dune::MPIHelper::getLocalCommunicator(); 102 103 #if PARMETIS_MAJOR_VERSION >= 4 104 const int OK = 105 #endif 106 ParMETIS_V3_PartMeshKway(elmdist.data(), eptr.data(), eind.data(), NULL, &wgtflag, &numflag, 107 &ncon, &ncommonnodes, &nparts, tpwgts.data(), ubvec.data(), 108 options, &edgecut, reinterpret_cast<idx_type*>(part.data()), &comm); 109 110 #if PARMETIS_MAJOR_VERSION >= 4 111 if (OK != METIS_OK) 112 DUNE_THROW(Dune::Exception, "ParMETIS returned an error code."); 113 #endif 114 } 115 116 return part; 117 } 118 119 /** \brief Create a repartitioning of a distributed Dune grid 120 * 121 * This pipes a Dune grid into the method ParMETIS_V3_AdaptiveRepart (see the ParMetis documentation for details) 122 * 123 * \param gv The grid view to be partitioned 124 * \param mpihelper The MPIHelper object, needed to get the MPI communicator 125 * \param itr ParMetis parameter to balance grid transport cost vs. communication cost for the redistributed grid. 126 * 127 * \return std::vector with one uint per All_Partition element. For each Interior_Partition element, the entry is the 128 * number of the partition the element is assigned to. 129 */ repartitionDune::ParMetisGridPartitioner130 static std::vector<unsigned> repartition(const GridView& gv, const Dune::MPIHelper& mpihelper, real_type itr = 1000) { 131 132 // Create global index map 133 GlobalIndexSet<GridView> globalIndex(gv,0); 134 135 int numElements = std::distance(gv.template begin<0, Interior_Partition>(), 136 gv.template end<0, Interior_Partition>()); 137 138 std::vector<unsigned> interiorPart(numElements); 139 140 // Setup parameters for ParMETIS 141 idx_type wgtflag = 0; // we don't use weights 142 idx_type numflag = 0; // we are using C-style arrays 143 idx_type ncon = 1; // number of balance constraints 144 idx_type options[4] = {0, 0, 0, 0}; // use default values for random seed, output and coupling 145 idx_type edgecut; // will store number of edges cut by partition 146 idx_type nparts = mpihelper.size(); // number of parts equals number of processes 147 std::vector<real_type> tpwgts(ncon*nparts, 1./nparts); // load per subdomain and weight (same load on every process) 148 std::vector<real_type> ubvec(ncon, 1.05); // weight tolerance (same weight tolerance for every weight there is) 149 150 MPI_Comm comm = Dune::MPIHelper::getCommunicator(); 151 152 // Make the number of interior elements of each processor available to all processors 153 std::vector<int> offset(gv.comm().size()); 154 std::fill(offset.begin(), offset.end(), 0); 155 156 gv.comm().template allgather<int>(&numElements, 1, offset.data()); 157 158 // The difference vtxdist[i+1] - vtxdist[i] is the number of elements that are on process i 159 std::vector<idx_type> vtxdist(gv.comm().size()+1); 160 vtxdist[0] = 0; 161 162 for (unsigned int i=1; i<vtxdist.size(); ++i) 163 vtxdist[i] = vtxdist[i-1] + offset[i-1]; 164 165 // Set up element adjacency lists 166 std::vector<idx_type> xadj, adjncy; 167 xadj.push_back(0); 168 169 for (const auto& element : elements(gv, Partitions::interior)) { 170 size_t numNeighbors = 0; 171 172 for (const auto& in : intersections(gv, element)) { 173 if (in.neighbor()) { 174 adjncy.push_back(globalIndex.index(in.outside())); 175 176 ++numNeighbors; 177 } 178 } 179 180 xadj.push_back(xadj.back() + numNeighbors); 181 } 182 183 #if PARMETIS_MAJOR_VERSION >= 4 184 const int OK = 185 #endif 186 ParMETIS_V3_AdaptiveRepart(vtxdist.data(), xadj.data(), adjncy.data(), NULL, NULL, NULL, 187 &wgtflag, &numflag, &ncon, &nparts, tpwgts.data(), ubvec.data(), 188 &itr, options, &edgecut, reinterpret_cast<idx_type*>(interiorPart.data()), &comm); 189 190 #if PARMETIS_MAJOR_VERSION >= 4 191 if (OK != METIS_OK) 192 DUNE_THROW(Dune::Exception, "ParMETIS returned error code " << OK); 193 #endif 194 195 // At this point, interiorPart contains a target rank for each interior element, and they are sorted 196 // by the order in which the grid view traverses them. Now we need to do two things: 197 // a) Add additional dummy entries for the ghost elements 198 // b) Use the element index for the actual ordering. Since there may be different types of elements, 199 // we cannot use the index set directly, but have to go through a Mapper. 200 201 typedef MultipleCodimMultipleGeomTypeMapper<GridView> ElementMapper; 202 ElementMapper elementMapper(gv, mcmgElementLayout()); 203 204 std::vector<unsigned int> part(gv.size(0)); 205 std::fill(part.begin(), part.end(), 0); 206 unsigned int c = 0; 207 for (const auto& element : elements(gv, Partitions::interior)) 208 part[elementMapper.index(element)] = interiorPart[c++]; 209 210 return part; 211 } 212 }; 213 214 } // namespace Dune 215 216 #else // PARMETIS_MAJOR_VERSION 217 #warning "You seem to be using the ParMETIS emulation layer of scotch, which does not work with this file." 218 #endif 219 220 #else // HAVE_PARMETIS 221 #warning "PARMETIS was not found, please check your configuration" 222 #endif 223 224 #endif // DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH 225