Line data Source code
1 : /*
2 : * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include "load_balancing.h"
34 : #include "lib_grid/algorithms/graph/dual_graph.h"
35 :
36 : // if we're using metis, then include the header now
37 : #ifdef UG_PARALLEL
38 : #include "pcl/pcl_base.h"
39 : #include "pcl/pcl_interface_communicator.h"
40 : #include "lib_grid/parallelization/util/parallel_dual_graph.h"
41 : #include "lib_grid/parallelization/util/compol_subset.h"
42 : #endif
43 :
44 : #ifdef UG_METIS
45 : extern "C" {
46 : #include "metis.h"
47 : }
48 : #endif
49 :
50 : #ifdef UG_PARMETIS
51 : extern "C" {
52 : #include "parmetis.h"
53 : }
54 : #endif
55 :
56 :
57 : using namespace std;
58 :
59 :
60 : namespace ug{
61 :
62 : ////////////////////////////////////////////////////////////////////////////////
63 : template <class TGeomBaseObj>
64 0 : bool PartitionGrid_MetisKway(SubsetHandler& shPartitionOut,
65 : Grid& grid, int numParts)
66 : {
67 : #ifdef UG_METIS
68 : if(numParts == 1){
69 : shPartitionOut.assign_subset(grid.begin<TGeomBaseObj>(),
70 : grid.end<TGeomBaseObj>(), 0);
71 : return true;
72 : }
73 :
74 : // construct the dual graph to the grid
75 : vector<idx_t> adjacencyMapStructure;
76 : vector<idx_t> adjacencyMap;
77 : ConstructDualGraph<TGeomBaseObj, idx_t>(adjacencyMapStructure,
78 : adjacencyMap, grid);
79 :
80 : // partition the graph using metis
81 : // first define options for the partitioning.
82 : idx_t options[METIS_NOPTIONS];
83 : METIS_SetDefaultOptions(options);
84 : options[METIS_OPTION_PTYPE] = METIS_PTYPE_KWAY;
85 : options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT;
86 : options[METIS_OPTION_NUMBERING] = 0;
87 : //request contiguous partitions
88 : //options[METIS_OPTION_CONTIG] = 1;
89 : //note: using the option METIS_OPTION_DBGLVL could be useful for debugging.
90 :
91 : int nVrts = (int)adjacencyMapStructure.size() - 1;
92 : int nConstraints = 1;
93 : int edgeCut;
94 : vector<idx_t> partitionMap(nVrts);
95 :
96 : UG_DLOG(LIB_GRID, 1, "CALLING METIS\n");
97 : int metisRet = METIS_PartGraphKway(&nVrts, &nConstraints,
98 : &adjacencyMapStructure.front(),
99 : &adjacencyMap.front(),
100 : NULL, NULL, NULL,
101 : &numParts, NULL, NULL, options,
102 : &edgeCut, &partitionMap.front());
103 : UG_DLOG(LIB_GRID, 1, "METIS DONE\n");
104 :
105 : if(metisRet != METIS_OK){
106 : UG_DLOG(LIB_GRID, 1, "METIS FAILED\n");
107 : return false;
108 : }
109 :
110 : // assign the subsets to the subset-handler
111 : int counter = 0;
112 : typedef typename geometry_traits<TGeomBaseObj>::iterator iterator;
113 : for(iterator iter = grid.begin<TGeomBaseObj>();
114 : iter != grid.end<TGeomBaseObj>(); ++iter)
115 : {
116 : shPartitionOut.assign_subset(*iter, partitionMap[counter++]);
117 : }
118 :
119 : return true;
120 : #else
121 0 : UG_THROW("METIS is not available in "
122 : "the current build. Please consider recompiling with METIS "
123 : "support enabled.\n");
124 : return false;
125 : #endif
126 : }
127 :
128 : //////////////////////////////
129 : // explicit instantiation
130 : template bool PartitionGrid_MetisKway<Edge>(SubsetHandler&, Grid&, int);
131 : template bool PartitionGrid_MetisKway<Face>(SubsetHandler&, Grid&, int);
132 : template bool PartitionGrid_MetisKway<Volume>(SubsetHandler&, Grid&, int);
133 :
134 :
135 : ////////////////////////////////////////////////////////////////////////////////
136 : template <class TGeomBaseObj>
137 0 : bool PartitionMultiGrid_MetisKway(SubsetHandler& shPartitionOut,
138 : MultiGrid& mg, int numParts,
139 : size_t baseLevel,
140 : int hWeight, int vWeight)
141 : {
142 : #ifdef UG_METIS
143 : // only call metis if more than 1 part is required
144 : if(numParts > 1){
145 : // here we'll store the dual graph
146 : vector<idx_t> adjacencyMapStructure;
147 : vector<idx_t> adjacencyMap;
148 : vector<idx_t> edgeWeightMap;
149 :
150 : //todo add baseLevel to ConstructDualGraphMG.
151 : ConstructDualGraphMG<TGeomBaseObj, idx_t>(adjacencyMapStructure,
152 : adjacencyMap, &edgeWeightMap,
153 : mg, baseLevel, hWeight, vWeight);
154 :
155 : // partition the graph using metis
156 : // first define options for the partitioning.
157 : idx_t options[METIS_NOPTIONS];
158 : METIS_SetDefaultOptions(options);
159 : options[METIS_OPTION_PTYPE] = METIS_PTYPE_KWAY;
160 : options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT;
161 : options[METIS_OPTION_NUMBERING] = 0;
162 : //request contiguous partitions
163 : //options[METIS_OPTION_CONTIG] = 1;
164 : //note: using the option METIS_OPTION_DBGLVL could be useful for debugging.
165 :
166 : int nVrts = (int)adjacencyMapStructure.size() - 1;
167 : int nConstraints = 1;
168 : int edgeCut;
169 : vector<idx_t> partitionMap(nVrts);
170 :
171 : UG_DLOG(LIB_GRID, 1, "CALLING METIS\n");
172 : int metisRet = METIS_PartGraphKway(&nVrts, &nConstraints,
173 : &adjacencyMapStructure.front(),
174 : &adjacencyMap.front(),
175 : NULL, NULL, &edgeWeightMap.front(),
176 : &numParts, NULL, NULL, options,
177 : &edgeCut, &partitionMap.front());
178 : UG_DLOG(LIB_GRID, 1, "METIS DONE\n");
179 :
180 : if(metisRet != METIS_OK){
181 : UG_DLOG(LIB_GRID, 1, "METIS FAILED\n");
182 : return false;
183 : }
184 :
185 :
186 : // assign the subsets to the subset-handler
187 : int counter = 0;
188 : for(size_t lvl = baseLevel; lvl < mg.num_levels(); ++lvl){
189 : typedef typename geometry_traits<TGeomBaseObj>::iterator iterator;
190 : for(iterator iter = mg.begin<TGeomBaseObj>(lvl);
191 : iter != mg.end<TGeomBaseObj>(lvl); ++iter)
192 : {
193 : shPartitionOut.assign_subset(*iter, partitionMap[counter++]);
194 : }
195 : }
196 : }
197 : else{
198 : // assign all elements to subset 0.
199 : for(size_t lvl = baseLevel; lvl < mg.num_levels(); ++lvl){
200 : shPartitionOut.assign_subset(mg.begin<TGeomBaseObj>(lvl),
201 : mg.end<TGeomBaseObj>(lvl), 0);
202 : }
203 : }
204 : return true;
205 : #else
206 0 : UG_THROW("METIS is not available in "
207 : "the current build. Please consider recompiling with METIS "
208 : "support enabled.\n");
209 : return false;
210 : #endif
211 : }
212 :
213 : ////////////////////////////////////////////////////////////////////////////////
214 : template <class TGeomBaseObj>
215 0 : bool PartitionMultiGrid_MetisKway(SubsetHandler& shPartitionOut,
216 : MultiGrid& mg, int numParts, size_t baseLevel,
217 : boost::function<int (TGeomBaseObj*, TGeomBaseObj*)>& weightFct)
218 : {
219 : #ifdef UG_METIS
220 : typedef TGeomBaseObj TElem;
221 : // only call metis if more than 1 part is required
222 : //if(numParts > 1){
223 : // here we'll store the dual graph
224 : vector<idx_t> adjacencyMapStructure;
225 : vector<idx_t> adjacencyMap;
226 : vector<idx_t> edgeWeightMap;
227 :
228 : if (weightFct)
229 : {
230 : // find out how many elems there are
231 : size_t size = 0;
232 : for (size_t lvl = baseLevel; lvl < mg.num_levels(); lvl++)
233 : size += mg.num<TElem>(lvl);
234 :
235 : // prepare list of all elements
236 : vector<TElem*> elems(size);
237 :
238 : // construct dual graph with standard edge weights
239 : ConstructDualGraphMG<TGeomBaseObj, idx_t>(adjacencyMapStructure,
240 : adjacencyMap, &edgeWeightMap,
241 : mg, baseLevel, 1, 1,
242 : NULL, &elems.front());
243 :
244 : // correct edge weights by calling weightFct for all pairs of TElems
245 : // found by ConstructDualGraphMG
246 : UG_ASSERT(size+1 == adjacencyMapStructure.size(), "index lengths not matching");
247 : TElem* elem1;
248 : TElem* elem2;
249 : for (size_t el = 0; el < size; el++)
250 : {
251 : elem1 = elems[el];
252 : for (size_t edgeInd = (size_t) adjacencyMapStructure[el]; edgeInd < (size_t) adjacencyMapStructure[el+1]; edgeInd++)
253 : {
254 : elem2 = elems[(size_t) adjacencyMap[edgeInd]];
255 : edgeWeightMap[edgeInd] = weightFct(elem1, elem2);
256 : }
257 : }
258 : }
259 : else
260 : {
261 : //todo add baseLevel to ConstructDualGraphMG.
262 : ConstructDualGraphMG<TGeomBaseObj, idx_t>(adjacencyMapStructure,
263 : adjacencyMap, &edgeWeightMap,
264 : mg, baseLevel, 1, 1);
265 : }
266 : if(numParts > 1){
267 : // partition the graph using metis
268 : // first define options for the partitioning.
269 : idx_t options[METIS_NOPTIONS];
270 : METIS_SetDefaultOptions(options);
271 : options[METIS_OPTION_PTYPE] = METIS_PTYPE_KWAY;
272 : options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT;
273 : options[METIS_OPTION_NUMBERING] = 0;
274 : //request contiguous partitions
275 : //options[METIS_OPTION_CONTIG] = 1;
276 : //note: using the option METIS_OPTION_DBGLVL could be useful for debugging.
277 :
278 : int nVrts = (int)adjacencyMapStructure.size() - 1;
279 : int nConstraints = 1;
280 : int edgeCut;
281 : vector<idx_t> partitionMap(nVrts);
282 :
283 : UG_DLOG(LIB_GRID, 1, "CALLING METIS\n");
284 : int metisRet = METIS_PartGraphKway(&nVrts, &nConstraints,
285 : &adjacencyMapStructure.front(),
286 : &adjacencyMap.front(),
287 : NULL, NULL, &edgeWeightMap.front(),
288 : &numParts, NULL, NULL, options,
289 : &edgeCut, &partitionMap.front());
290 : UG_DLOG(LIB_GRID, 1, "METIS DONE\n");
291 :
292 : if(metisRet != METIS_OK){
293 : UG_DLOG(LIB_GRID, 1, "METIS FAILED\n");
294 : return false;
295 : }
296 :
297 :
298 : // assign the subsets to the subset-handler
299 : int counter = 0;
300 : for(size_t lvl = baseLevel; lvl < mg.num_levels(); ++lvl){
301 : typedef typename geometry_traits<TGeomBaseObj>::iterator iterator;
302 : for(iterator iter = mg.begin<TGeomBaseObj>(lvl);
303 : iter != mg.end<TGeomBaseObj>(lvl); ++iter)
304 : {
305 : shPartitionOut.assign_subset(*iter, partitionMap[counter++]);
306 : }
307 : }
308 : }
309 : else{
310 : // assign all elements to subset 0.
311 : for(size_t lvl = baseLevel; lvl < mg.num_levels(); ++lvl){
312 : shPartitionOut.assign_subset(mg.begin<TGeomBaseObj>(lvl),
313 : mg.end<TGeomBaseObj>(lvl), 0);
314 : }
315 : }
316 : return true;
317 : #else
318 0 : UG_THROW("METIS is not available in "
319 : "the current build. Please consider recompiling with METIS "
320 : "support enabled.\n");
321 : return false;
322 : #endif
323 : }
324 :
325 : ////////////////////////////////////////////////////////////////////////////////
326 : template <class TGeomBaseObj>
327 0 : bool PartitionMultiGridLevel_MetisKway(SubsetHandler& shPartitionOut,
328 : MultiGrid& mg, int numParts, size_t level)
329 : {
330 0 : if(level > mg.top_level()){
331 : UG_LOG(" WARNING in PartitionMultiGridLevel_MetisKway:"
332 : "Specified level too high: toplevel=" << mg.top_level()
333 : << ", specified-level=" << level <<"\n");
334 0 : return true;
335 : }
336 :
337 0 : if(mg.num<TGeomBaseObj>() == 0)
338 : return true;
339 :
340 : #ifdef UG_METIS
341 : typedef TGeomBaseObj TElem;
342 : typedef typename geometry_traits<TGeomBaseObj>::iterator ElemIter;
343 :
344 : // only call metis if more than 1 part is required
345 : int rootProc = 0;
346 : #ifdef UG_PARALLEL
347 : rootProc = pcl::ProcRank();
348 : #endif
349 :
350 : if(numParts > 1){
351 : // here we'll store the dual graph
352 : vector<idx_t> adjacencyMapStructure;
353 : vector<idx_t> adjacencyMap;
354 :
355 :
356 : // we can optionally use a higher edge weight on siblings.
357 : // note - higher memory and processing requirements if enabled.
358 : const bool useEdgeWeights = true;
359 : const idx_t siblingWeight = 2;
360 :
361 : // we'll reuse the index later on during assignment of the edge weights
362 : typedef Attachment<idx_t> AIndex;
363 : AIndex aIndex;
364 : vector<TElem*> elems;
365 :
366 : if(useEdgeWeights){
367 : mg.attach_to<TElem>(aIndex);
368 : Grid::AttachmentAccessor<TElem, AIndex> aaInd;
369 : elems.resize(mg.num<TElem>(level));
370 :
371 : // we also need a vector of elements, so that we can access them by index
372 : // (also for the caclulation of edge weights)
373 : ConstructDualGraphMGLevel<TElem, idx_t>(adjacencyMapStructure, adjacencyMap,
374 : mg, level, &aIndex, &elems.front());
375 : }
376 : else{
377 : ConstructDualGraphMGLevel<TElem, idx_t>(adjacencyMapStructure, adjacencyMap,
378 : mg, level);
379 : }
380 :
381 : // partition the graph using metis
382 : // first define options for the partitioning.
383 : idx_t options[METIS_NOPTIONS];
384 : METIS_SetDefaultOptions(options);
385 : options[METIS_OPTION_PTYPE] = METIS_PTYPE_KWAY;
386 : options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT;
387 : options[METIS_OPTION_NUMBERING] = 0;
388 : //request contiguous partitions
389 : //options[METIS_OPTION_CONTIG] = 1;
390 : //note: using the option METIS_OPTION_DBGLVL could be useful for debugging.
391 :
392 : idx_t nVrts = (idx_t)adjacencyMapStructure.size() - 1;
393 : idx_t nConstraints = 1;
394 : idx_t edgeCut;
395 : vector<idx_t> partitionMap(nVrts);
396 :
397 : // create a weight map for the vertices based on the number of children+1
398 : // for each graph-vertex. This is not necessary, if we're already on the top level
399 : idx_t* pVrtSizeMap = NULL;
400 : vector<idx_t> vrtSizeMap;
401 : if(level < mg.top_level()){
402 : vrtSizeMap.reserve(nVrts);
403 : for(ElemIter iter = mg.begin<TElem>(level);
404 : iter != mg.end<TElem>(level); ++iter)
405 : {
406 : vrtSizeMap.push_back((mg.num_children_total(*iter) + 1));
407 : }
408 : assert((int)vrtSizeMap.size() == nVrts);
409 : pVrtSizeMap = &vrtSizeMap.front();
410 : }
411 :
412 : // we'll also create weights for the edges, since we want to cluster elements,
413 : // which have the same parent (this reduces several problems later on, like
414 : // additional vertical interfaces)
415 : idx_t* pEdgeWeights = NULL;
416 : vector<idx_t> edgeWeights;
417 : if(useEdgeWeights){
418 : edgeWeights.reserve(adjacencyMap.size());
419 : for(int i_vrt = 0; i_vrt < nVrts; ++i_vrt){
420 : int start = adjacencyMapStructure[i_vrt];
421 : int end = adjacencyMapStructure[i_vrt + 1];
422 :
423 : GridObject* parent = mg.get_parent(elems[i_vrt]);
424 : for(int i_edge = start; i_edge < end; ++i_edge){
425 : if(parent == mg.get_parent(elems[adjacencyMap[i_edge]]))
426 : edgeWeights.push_back(siblingWeight);
427 : else
428 : edgeWeights.push_back(1);
429 : }
430 : }
431 : pEdgeWeights = &edgeWeights.front();
432 : }
433 :
434 : UG_DLOG(LIB_GRID, 1, "CALLING METIS\n");
435 : int metisRet = METIS_PartGraphKway(&nVrts, &nConstraints,
436 : &adjacencyMapStructure.front(),
437 : &adjacencyMap.front(),
438 : NULL, pVrtSizeMap, pEdgeWeights,
439 : &numParts, NULL, NULL, options,
440 : &edgeCut, &partitionMap.front());
441 : UG_DLOG(LIB_GRID, 1, "METIS DONE\n");
442 :
443 : if(metisRet != METIS_OK){
444 : UG_DLOG(LIB_GRID, 1, "METIS FAILED\n");
445 : return false;
446 : }
447 :
448 :
449 : // assign the subsets to the subset-handler
450 : // all subsets below the specified level go to the rootProc.
451 : // The ones on the specified level go as METIS assigned them.
452 : // All children in levels above copy their from their parent.
453 : for(size_t lvl = 0; lvl < level; ++lvl)
454 : for(ElemIter iter = mg.begin<TElem>(lvl); iter != mg.end<TElem>(lvl); ++iter)
455 : shPartitionOut.assign_subset(*iter, rootProc);
456 :
457 : int counter = 0;
458 : for(ElemIter iter = mg.begin<TElem>(level); iter != mg.end<TElem>(level); ++iter)
459 : shPartitionOut.assign_subset(*iter, partitionMap[counter++]);
460 :
461 : // currently there is a restriction in the functionality of the surface view.
462 : // Because of that we have to make sure, that all siblings in the specified level
463 : // are sent to the same process... we thus adjust the partition slightly.
464 : //todo: Not all siblings should have to be sent to the same process...
465 : // simply remove the following code block - make sure that surface-view supports this!
466 : // However, problems with discretizations and solvers would occur
467 : if(level > 0){
468 : // put all children in the subset of the first one.
469 : for(ElemIter iter = mg.begin<TElem>(level-1);
470 : iter != mg.end<TElem>(level-1); ++iter)
471 : {
472 : TElem* e = *iter;
473 : size_t numChildren = mg.num_children<TElem>(e);
474 : if(numChildren > 1){
475 : int partition = shPartitionOut.get_subset_index(mg.get_child<TElem>(e, 0));
476 : for(size_t i = 1; i < numChildren; ++i)
477 : shPartitionOut.assign_subset(mg.get_child<TElem>(e, i), partition);
478 : }
479 : }
480 : }
481 :
482 : for(size_t lvl = level; lvl < mg.top_level(); ++lvl){
483 : for(ElemIter iter = mg.begin<TElem>(lvl); iter != mg.end<TElem>(lvl); ++iter)
484 : {
485 : size_t numChildren = mg.num_children<TElem>(*iter);
486 : for(size_t i = 0; i < numChildren; ++i){
487 : shPartitionOut.assign_subset(mg.get_child<TElem>(*iter, i),
488 : shPartitionOut.get_subset_index(*iter));
489 : }
490 : }
491 : }
492 :
493 : if(useEdgeWeights){
494 : mg.detach_from<TElem>(aIndex);
495 : }
496 : }
497 : else{
498 : // assign all elements to subset 0.
499 : for(size_t lvl = 0; lvl < mg.num_levels(); ++lvl){
500 : shPartitionOut.assign_subset(mg.begin<TGeomBaseObj>(lvl),
501 : mg.end<TGeomBaseObj>(lvl), rootProc);
502 : }
503 : }
504 : return true;
505 : #else
506 0 : UG_THROW("METIS is not available in "
507 : "the current build. Please consider recompiling with METIS "
508 : "support enabled.\n");
509 : return false;
510 : #endif
511 : }
512 :
513 : ////////////////////////////////////////////////////////////////////////////////
514 : template <class TGeomBaseObj>
515 0 : bool PartitionMultiGridLevel_ParmetisKway(SubsetHandler& shPartitionOut,
516 : MultiGrid& mg, int numParts, size_t level)
517 : {
518 : UG_DLOG(LIB_GRID, 1, "start - PartitionMultiGridLevel_ParmetisKway\n");
519 :
520 : #if defined UG_PARMETIS && defined UG_PARALLEL
521 : typedef TGeomBaseObj TElem;
522 : typedef typename geometry_traits<TGeomBaseObj>::iterator ElemIter;
523 :
524 : int localProc = pcl::ProcRank();
525 : pcl::ProcessCommunicator procComWorld;
526 :
527 : // here we'll store the dual graph
528 : vector<idx_t> adjacencyMapStructure;
529 : vector<idx_t> adjacencyMap;
530 : vector<idx_t> nodeOffsetMap;
531 :
532 : ConstructParallelDualGraphMGLevel<TElem, idx_t>(adjacencyMapStructure,
533 : adjacencyMap, nodeOffsetMap,
534 : mg, level, procComWorld);
535 :
536 : UG_DLOG(LIB_GRID, 2, " parallel dual graph #vrts: " << (int)adjacencyMapStructure.size() - 1
537 : << ", #edges: " << adjacencyMap.size() / 2 << "\n");
538 :
539 : pcl::ProcessCommunicator procCom = procComWorld.create_sub_communicator(adjacencyMap.size() > 1);
540 :
541 : /*if(!pcl::AllProcsTrue(adjacencyMap.size() > 1)){
542 : UG_THROW("ParMetis may only be executed if all processes contain non-ghost elements"
543 : " on the given level (at least two neighboring).")
544 : }*/
545 :
546 : // partition the graph using parmetis
547 : idx_t options[3]; options[0] = 0;//default values
548 : idx_t nVrts = (idx_t)adjacencyMapStructure.size() - 1;
549 : idx_t nConstraints = 1;
550 : idx_t edgeCut;
551 : idx_t wgtFlag = 2;//only vertices are weighted
552 : idx_t numFlag = 0;
553 : vector<idx_t> partitionMap(nVrts);
554 : vector<real_t> tpwgts(numParts, 1. / (number)numParts);
555 : real_t ubvec = 1.05;
556 : // create a weight map for the vertices based on the number of children+1
557 : // for each graph-vertex. This is not necessary, if we're already on the top level
558 : idx_t* pVrtSizeMap = NULL;
559 : vector<idx_t> vrtSizeMap;
560 : {
561 : vrtSizeMap.reserve(nVrts);
562 : for(ElemIter iter = mg.begin<TElem>(level);
563 : iter != mg.end<TElem>(level); ++iter)
564 : {
565 : vrtSizeMap.push_back((mg.num_children_total(*iter) + 1));
566 : }
567 : assert((int)vrtSizeMap.size() == nVrts);
568 : pVrtSizeMap = &vrtSizeMap.front();
569 : }
570 :
571 : if(!procCom.empty()){
572 : UG_DLOG(LIB_GRID, 1, "CALLING PARMETIS\n");
573 : MPI_Comm mpiCom = procCom.get_mpi_communicator();
574 : int metisRet = ParMETIS_V3_PartKway(&nodeOffsetMap.front(),
575 : &adjacencyMapStructure.front(),
576 : &adjacencyMap.front(),
577 : pVrtSizeMap, NULL, &wgtFlag,
578 : &numFlag, &nConstraints,
579 : &numParts, &tpwgts.front(), &ubvec, options,
580 : &edgeCut, &partitionMap.front(),
581 : &mpiCom);
582 : UG_DLOG(LIB_GRID, 1, "PARMETIS DONE\n");
583 :
584 : if(metisRet != METIS_OK){
585 : UG_THROW("PARMETIS FAILED on process " << localProc);
586 : }
587 : }
588 :
589 : // assign the subsets to the subset-handler
590 : // all subsets below the specified level go to the rootProc.
591 : // The ones on the specified level go as METIS assigned them.
592 : // All children in levels above copy their from their parent.
593 : for(size_t lvl = 0; lvl < level; ++lvl)
594 : for(ElemIter iter = mg.begin<TElem>(lvl); iter != mg.end<TElem>(lvl); ++iter)
595 : shPartitionOut.assign_subset(*iter, localProc);
596 :
597 : assert(mg.is_parallel());
598 : DistributedGridManager& distGridMgr = *mg.distributed_grid_manager();
599 : int counter = 0;
600 : for(ElemIter iter = mg.begin<TElem>(level); iter != mg.end<TElem>(level); ++iter){
601 : if(!distGridMgr.is_ghost(*iter))
602 : shPartitionOut.assign_subset(*iter, partitionMap[counter++]);
603 : }
604 :
605 :
606 :
607 : typedef typename GridLayoutMap::Types<TElem>::Layout::LevelLayout ElemLayout;
608 : GridLayoutMap& glm = mg.distributed_grid_manager()->grid_layout_map();
609 : pcl::InterfaceCommunicator<ElemLayout> com;
610 : ComPol_Subset<ElemLayout> compolSHCopy(shPartitionOut, true);
611 :
612 : // copy subset indices from vertical slaves to vertical masters,
613 : // since partitioning was only performed on vslaves
614 : if(glm.has_layout<TElem>(INT_V_SLAVE))
615 : com.send_data(glm.get_layout<TElem>(INT_V_SLAVE).layout_on_level(level),
616 : compolSHCopy);
617 : if(glm.has_layout<TElem>(INT_V_MASTER))
618 : com.receive_data(glm.get_layout<TElem>(INT_V_MASTER).layout_on_level(level),
619 : compolSHCopy);
620 : com.communicate();
621 :
622 : //todo: make sure that there are no problems at vertical interfaces
623 : // currently there is a restriction in the functionality of the surface view.
624 : // Because of that we have to make sure, that all siblings in the specified level
625 : // are sent to the same process... we thus adjust the partition slightly.
626 : //todo: Not all siblings should have to be sent to the same process...
627 : // simply remove the following code block - make sure that surface-view supports this!
628 : // However, problems with discretizations and solvers would occur.
629 : // If you remove the following block, consider reducing v-communication.
630 : if(level > 0){
631 : // put all children in the subset of the first one.
632 : for(ElemIter iter = mg.begin<TElem>(level-1);
633 : iter != mg.end<TElem>(level-1); ++iter)
634 : {
635 : TElem* e = *iter;
636 : size_t numChildren = mg.num_children<TElem>(e);
637 : if(numChildren > 1){
638 : int partition = shPartitionOut.get_subset_index(mg.get_child<TElem>(e, 0));
639 : for(size_t i = 1; i < numChildren; ++i)
640 : shPartitionOut.assign_subset(mg.get_child<TElem>(e, i), partition);
641 : }
642 : }
643 : }
644 :
645 :
646 : for(size_t lvl = level; lvl < mg.top_level(); ++lvl){
647 : // copy subset indices from vertical masters to vertical slaves
648 : if(glm.has_layout<TElem>(INT_V_MASTER))
649 : com.send_data(glm.get_layout<TElem>(INT_V_MASTER).layout_on_level(lvl),
650 : compolSHCopy);
651 : if(glm.has_layout<TElem>(INT_V_SLAVE))
652 : com.receive_data(glm.get_layout<TElem>(INT_V_SLAVE).layout_on_level(lvl),
653 : compolSHCopy);
654 : com.communicate();
655 :
656 : for(ElemIter iter = mg.begin<TElem>(lvl); iter != mg.end<TElem>(lvl); ++iter)
657 : {
658 : size_t numChildren = mg.num_children<TElem>(*iter);
659 : for(size_t i = 0; i < numChildren; ++i){
660 : shPartitionOut.assign_subset(mg.get_child<TElem>(*iter, i),
661 : shPartitionOut.get_subset_index(*iter));
662 : }
663 : }
664 : }
665 :
666 : // and a final copy on the top level...
667 : if(mg.num_levels() > 1){
668 : if(glm.has_layout<TElem>(INT_V_MASTER))
669 : com.send_data(glm.get_layout<TElem>(INT_V_MASTER).layout_on_level(mg.top_level()),
670 : compolSHCopy);
671 : if(glm.has_layout<TElem>(INT_V_SLAVE))
672 : com.receive_data(glm.get_layout<TElem>(INT_V_SLAVE).layout_on_level(mg.top_level()),
673 : compolSHCopy);
674 : com.communicate();
675 : }
676 :
677 : UG_DLOG(LIB_GRID, 1, "stop - PartitionMultiGridLevel_ParmetisKway\n");
678 : return true;
679 : #else
680 0 : UG_THROW("METIS is not available in "
681 : "the current build. Please consider recompiling with METIS "
682 : "support enabled.\n");
683 : return false;
684 : #endif
685 : }
686 :
687 :
688 : //////////////////////////////
689 : // explicit instantiation
690 : template bool PartitionMultiGrid_MetisKway<Edge>(SubsetHandler&, MultiGrid&,
691 : int, size_t, int, int);
692 : template bool PartitionMultiGrid_MetisKway<Face>(SubsetHandler&, MultiGrid&,
693 : int, size_t, int, int);
694 : template bool PartitionMultiGrid_MetisKway<Volume>(SubsetHandler&, MultiGrid&,
695 : int, size_t, int, int);
696 :
697 : template bool PartitionMultiGrid_MetisKway<Edge>(SubsetHandler&, MultiGrid&, int, size_t,
698 : boost::function<int (Edge*, Edge*)>&);
699 : template bool PartitionMultiGrid_MetisKway<Face>(SubsetHandler&, MultiGrid&, int, size_t,
700 : boost::function<int (Face*, Face*)>&);
701 : template bool PartitionMultiGrid_MetisKway<Volume>(SubsetHandler&, MultiGrid&, int, size_t,
702 : boost::function<int (Volume*, Volume*)>&);
703 :
704 : template bool PartitionMultiGridLevel_MetisKway<Edge>(SubsetHandler&,
705 : MultiGrid&, int, size_t);
706 : template bool PartitionMultiGridLevel_MetisKway<Face>(SubsetHandler&,
707 : MultiGrid&, int, size_t);
708 : template bool PartitionMultiGridLevel_MetisKway<Volume>(SubsetHandler&,
709 : MultiGrid&, int, size_t);
710 :
711 : template bool PartitionMultiGridLevel_ParmetisKway<Edge>(SubsetHandler&,
712 : MultiGrid&, int, size_t);
713 : template bool PartitionMultiGridLevel_ParmetisKway<Face>(SubsetHandler&,
714 : MultiGrid&, int, size_t);
715 : template bool PartitionMultiGridLevel_ParmetisKway<Volume>(SubsetHandler&,
716 : MultiGrid&, int, size_t);
717 :
718 : }// end of namespace
|
