Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local Includes
21 : #include "libmesh/metis_partitioner.h"
22 :
23 : #include "libmesh/libmesh_config.h"
24 : #include "libmesh/elem.h"
25 : #include "libmesh/enum_partitioner_type.h"
26 : #include "libmesh/error_vector.h"
27 : #include "libmesh/libmesh_logging.h"
28 : #include "libmesh/mesh_base.h"
29 : #include "libmesh/mesh_communication.h"
30 : #include "libmesh/mesh_tools.h"
31 : #include "libmesh/metis_csr_graph.h"
32 : #include "libmesh/parallel.h"
33 : #include "libmesh/utility.h"
34 : #include "libmesh/vectormap.h"
35 :
36 : #ifdef LIBMESH_HAVE_METIS
37 : // MIPSPro 7.4.2 gets confused about these nested namespaces
38 : # ifdef __sgi
39 : # include <cstdarg>
40 : # endif
41 : namespace Metis {
42 : extern "C" {
43 : # include "libmesh/ignore_warnings.h"
44 : # include "metis.h"
45 : # include "libmesh/restore_warnings.h"
46 : }
47 : }
48 : #else
49 : # include "libmesh/sfc_partitioner.h"
50 : #endif
51 :
52 :
53 : // C++ includes
54 : #include <map>
55 : #include <unordered_map>
56 : #include <vector>
57 :
58 :
59 : namespace libMesh
60 : {
61 :
62 :
63 142 : PartitionerType MetisPartitioner::type() const
64 : {
65 142 : return METIS_PARTITIONER;
66 : }
67 :
68 :
69 233102 : void MetisPartitioner::partition_range(MeshBase & mesh,
70 : MeshBase::element_iterator beg,
71 : MeshBase::element_iterator end,
72 : unsigned int n_pieces)
73 : {
74 : // Check for easy returns
75 233102 : if (beg == end)
76 71 : return;
77 :
78 233031 : if (n_pieces == 1)
79 : {
80 0 : this->single_partition_range (beg, end);
81 0 : return;
82 : }
83 :
84 8000 : libmesh_assert_greater (n_pieces, 0);
85 :
86 : // We don't yet support distributed meshes with this Partitioner
87 233031 : if (!mesh.is_serial())
88 : {
89 0 : libMesh::out << "WARNING: Forced to gather a distributed mesh for METIS" << std::endl;
90 0 : mesh.allgather();
91 : }
92 :
93 : // What to do if the Metis library IS NOT present
94 : #ifndef LIBMESH_HAVE_METIS
95 :
96 : libmesh_do_once(
97 : libMesh::out << "ERROR: The library has been built without" << std::endl
98 : << "Metis support. Using a space-filling curve" << std::endl
99 : << "partitioner instead!" << std::endl;);
100 :
101 : SFCPartitioner sfcp;
102 : sfcp.partition_range (mesh, beg, end, n_pieces);
103 :
104 : // What to do if the Metis library IS present
105 : #else
106 :
107 16000 : LOG_SCOPE("partition_range()", "MetisPartitioner");
108 :
109 458062 : const dof_id_type n_range_elem = cast_int<dof_id_type>(std::distance(beg, end));
110 :
111 : // Metis will only consider the elements in the range.
112 : // We need to map the range element ids into a
113 : // contiguous range. Further, we want the unique range indexing to be
114 : // independent of the element ordering, otherwise a circular dependency
115 : // can result in which the partitioning depends on the ordering which
116 : // depends on the partitioning...
117 16000 : vectormap<dof_id_type, dof_id_type> global_index_map;
118 233031 : global_index_map.reserve (n_range_elem);
119 :
120 : {
121 16000 : std::vector<dof_id_type> global_index;
122 :
123 233031 : MeshCommunication().find_global_indices (mesh.comm(),
124 233031 : MeshTools::create_bounding_box(mesh),
125 : beg, end, global_index);
126 :
127 8000 : libmesh_assert_equal_to (global_index.size(), n_range_elem);
128 :
129 : // If we have unique_id() then global_index entries should be
130 : // unique, so we just need a map for them.
131 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
132 8000 : std::size_t cnt=0;
133 13322829 : for (const auto & elem : as_range(beg, end))
134 24900037 : global_index_map.emplace(elem->id(), global_index[cnt++]);
135 : #else
136 : // If we don't have unique_id() then Hilbert SFC based
137 : // global_index entries might overlap if elements overlap, so we
138 : // need to fix any duplicates ...
139 : //
140 : // First, check for duplicates in O(N) time
141 : std::unordered_map<dof_id_type, std::vector<dof_id_type>>
142 : global_index_inverses;
143 : bool found_duplicate_indices = false;
144 : {
145 : std::size_t cnt=0;
146 : for (const auto & elem : as_range(beg, end))
147 : {
148 : auto & vec = global_index_inverses[global_index[cnt++]];
149 : if (!vec.empty())
150 : found_duplicate_indices = true;
151 : vec.push_back(elem->id());
152 : }
153 : }
154 :
155 : // But if we find them, we'll need O(N log N) to fix them while
156 : // retaining the same space-filling-curve for efficient initial
157 : // partitioning.
158 : if (found_duplicate_indices)
159 : {
160 : const std::map<dof_id_type, std::vector<dof_id_type>>
161 : sorted_inverses(global_index_inverses.begin(),
162 : global_index_inverses.end());
163 : std::size_t new_global_index=0;
164 : for (const auto & [index, elem_ids] : sorted_inverses)
165 : for (const dof_id_type elem_id : elem_ids)
166 : global_index_map.emplace(elem_id, new_global_index++);
167 : }
168 : else
169 : {
170 : std::size_t cnt=0;
171 : for (const auto & elem : as_range(beg, end))
172 : global_index_map.emplace(elem->id(), global_index[cnt++]);
173 : }
174 : #endif
175 :
176 8000 : libmesh_assert_equal_to (global_index_map.size(), n_range_elem);
177 :
178 : #ifdef DEBUG
179 16000 : std::unordered_set<dof_id_type> distinct_global_indices, distinct_ids;
180 1054528 : for (const auto & [elem_id, index] : global_index_map)
181 : {
182 1046528 : distinct_ids.insert(elem_id);
183 1046528 : distinct_global_indices.insert(index);
184 : }
185 8000 : libmesh_assert_equal_to (distinct_global_indices.size(), n_range_elem);
186 8000 : libmesh_assert_equal_to (distinct_ids.size(), n_range_elem);
187 : #endif
188 : }
189 :
190 : // If we have boundary elements in this mesh, we want to account for
191 : // the connectivity between them and interior elements. We can find
192 : // interior elements from boundary elements, but we need to build up
193 : // a lookup map to do the reverse.
194 : typedef std::unordered_multimap<const Elem *, const Elem *> map_type;
195 16000 : map_type interior_to_boundary_map;
196 :
197 24094538 : for (const auto & elem : as_range(beg, end))
198 : {
199 : // If we don't have an interior_parent then there's nothing
200 : // to look us up.
201 20645077 : if ((elem->dim() >= LIBMESH_DIM) ||
202 7780310 : !elem->interior_parent())
203 12861494 : continue;
204 :
205 : // get all relevant interior elements
206 880 : std::set<Elem *> neighbor_set;
207 3273 : elem->find_interior_neighbors(neighbor_set);
208 :
209 8685 : for (auto & neighbor : neighbor_set)
210 502 : interior_to_boundary_map.emplace(neighbor, elem);
211 217031 : }
212 :
213 : // Data structure that Metis will fill up on processor 0 and broadcast.
214 241031 : std::vector<Metis::idx_t> part(n_range_elem);
215 :
216 : // Invoke METIS, but only on processor 0.
217 : // Then broadcast the resulting decomposition
218 241031 : if (mesh.processor_id() == 0)
219 : {
220 : // Data structures and parameters needed only on processor 0 by Metis.
221 : // std::vector<Metis::idx_t> options(5);
222 37961 : std::vector<Metis::idx_t> vwgt(n_range_elem);
223 :
224 : Metis::idx_t
225 33961 : n = static_cast<Metis::idx_t>(n_range_elem), // number of "nodes" (elements) in the graph
226 : // wgtflag = 2, // weights on vertices only, none on edges
227 : // numflag = 0, // C-style 0-based numbering
228 33961 : nparts = static_cast<Metis::idx_t>(n_pieces), // number of subdomains to create
229 33961 : edgecut = 0; // the numbers of edges cut by the resulting partition
230 :
231 : // Set the options
232 : // options[0] = 0; // use default options
233 :
234 : // build the graph
235 12000 : METIS_CSR_Graph<Metis::idx_t> csr_graph;
236 :
237 33961 : csr_graph.offsets.resize(n_range_elem + 1, 0);
238 :
239 : // Local scope for these
240 : {
241 : // build the graph in CSR format. Note that
242 : // the edges in the graph will correspond to
243 : // face neighbors
244 :
245 : #ifdef LIBMESH_ENABLE_AMR
246 8000 : std::vector<const Elem *> neighbors_offspring;
247 : #endif
248 :
249 : #ifndef NDEBUG
250 4000 : std::size_t graph_size=0;
251 : #endif
252 :
253 : // (1) first pass - get the row sizes for each element by counting the number
254 : // of face neighbors. Also populate the vwght array if necessary
255 4101673 : for (const auto & elem : as_range(beg, end))
256 : {
257 : const dof_id_type elem_global_index =
258 2542140 : global_index_map[elem->id()];
259 :
260 523264 : libmesh_assert_less (elem_global_index, vwgt.size());
261 :
262 : // The weight is used to define what a balanced graph is
263 2542140 : if (!_weights)
264 : {
265 : // Spline nodes are a special case (storing all the
266 : // unconstrained DoFs in an IGA simulation), but in
267 : // general we'll try to distribute work by expecting
268 : // it to be roughly proportional to DoFs, which are
269 : // roughly proportional to nodes.
270 2544127 : if (elem->type() == NODEELEM &&
271 6667 : elem->mapping_type() == RATIONAL_BERNSTEIN_MAP)
272 7813 : vwgt[elem_global_index] = 50;
273 : else
274 2536154 : vwgt[elem_global_index] = elem->n_nodes();
275 : }
276 : else
277 0 : vwgt[elem_global_index] = static_cast<Metis::idx_t>((*_weights)[elem->id()]);
278 :
279 523264 : unsigned int num_neighbors = 0;
280 :
281 : // Loop over the element's neighbors. An element
282 : // adjacency corresponds to a face neighbor
283 13167127 : for (auto neighbor : elem->neighbor_ptr_range())
284 : {
285 10101722 : if (neighbor != nullptr)
286 : {
287 : // If the neighbor is active, but is not in the
288 : // range of elements being partitioned, treat it
289 : // as a nullptr neighbor.
290 9282750 : if (neighbor->active() && !global_index_map.count(neighbor->id()))
291 0 : continue;
292 :
293 : // If the neighbor is active treat it
294 : // as a connection
295 2005101 : if (neighbor->active())
296 9183949 : num_neighbors++;
297 :
298 : #ifdef LIBMESH_ENABLE_AMR
299 :
300 : // Otherwise we need to find all of the
301 : // neighbor's children that are connected to
302 : // us and add them
303 : else
304 : {
305 : // The side of the neighbor to which
306 : // we are connected
307 : const unsigned int ns =
308 98801 : neighbor->which_neighbor_am_i (elem);
309 18141 : libmesh_assert_less (ns, neighbor->n_neighbors());
310 :
311 : // Get all the active children (& grandchildren, etc...)
312 : // of the neighbor.
313 :
314 : // FIXME - this is the wrong thing, since we
315 : // should be getting the active family tree on
316 : // our side only. But adding too many graph
317 : // links may cause hanging nodes to tend to be
318 : // on partition interiors, which would reduce
319 : // communication overhead for constraint
320 : // equations, so we'll leave it.
321 98801 : neighbor->active_family_tree (neighbors_offspring);
322 :
323 : // Get all the neighbor's children that
324 : // live on that side and are thus connected
325 : // to us
326 584297 : for (const auto & child : neighbors_offspring)
327 : {
328 : // Skip neighbor offspring which are not in the range of elements being partitioned.
329 485496 : if (!global_index_map.count(child->id()))
330 0 : continue;
331 :
332 : // This does not assume a level-1 mesh.
333 : // Note that since children have sides numbered
334 : // coincident with the parent then this is a sufficient test.
335 579031 : if (child->neighbor_ptr(ns) == elem)
336 : {
337 40870 : libmesh_assert (child->active());
338 211079 : num_neighbors++;
339 : }
340 : }
341 : }
342 :
343 : #endif /* ifdef LIBMESH_ENABLE_AMR */
344 :
345 : }
346 : }
347 :
348 : // Check for any interior neighbors
349 2542140 : if ((elem->dim() < LIBMESH_DIM) && elem->interior_parent())
350 : {
351 : // get all relevant interior elements
352 220 : std::set<const Elem *> neighbor_set;
353 877 : elem->find_interior_neighbors(neighbor_set);
354 :
355 877 : num_neighbors += cast_int<unsigned int>(neighbor_set.size());
356 : }
357 :
358 : // Check for any boundary neighbors
359 : typedef map_type::iterator map_it_type;
360 : std::pair<map_it_type, map_it_type>
361 1046529 : bounds = interior_to_boundary_map.equal_range(elem);
362 2542140 : num_neighbors += cast_int<unsigned int>
363 523264 : (std::distance(bounds.first, bounds.second));
364 :
365 : // Whether we enabled unique_id or had to fix any overlaps
366 : // by hand, our global indices should be unique, so we
367 : // should never see the same elem_global_index twice, so
368 : // we should never be trying to edit an already-nonzero
369 : // offset.
370 523264 : libmesh_assert(!csr_graph.offsets[elem_global_index+1]);
371 :
372 523264 : csr_graph.prep_n_nonzeros(elem_global_index, num_neighbors);
373 : #ifndef NDEBUG
374 523264 : graph_size += num_neighbors;
375 : #endif
376 25961 : }
377 :
378 33961 : csr_graph.prepare_for_use();
379 :
380 : // (2) second pass - fill the compressed adjacency array
381 4101673 : for (const auto & elem : as_range(beg, end))
382 : {
383 : const dof_id_type elem_global_index =
384 2542140 : global_index_map[elem->id()];
385 :
386 523264 : unsigned int connection=0;
387 :
388 : // Loop over the element's neighbors. An element
389 : // adjacency corresponds to a face neighbor
390 13167127 : for (auto neighbor : elem->neighbor_ptr_range())
391 : {
392 10101722 : if (neighbor != nullptr)
393 : {
394 : // If the neighbor is active, but is not in the
395 : // range of elements being partitioned, treat it
396 : // as a nullptr neighbor.
397 9282750 : if (neighbor->active() && !global_index_map.count(neighbor->id()))
398 0 : continue;
399 :
400 : // If the neighbor is active treat it
401 : // as a connection
402 2005101 : if (neighbor->active())
403 9183949 : csr_graph(elem_global_index, connection++) = global_index_map[neighbor->id()];
404 :
405 : #ifdef LIBMESH_ENABLE_AMR
406 :
407 : // Otherwise we need to find all of the
408 : // neighbor's children that are connected to
409 : // us and add them
410 : else
411 : {
412 : // The side of the neighbor to which
413 : // we are connected
414 : const unsigned int ns =
415 98801 : neighbor->which_neighbor_am_i (elem);
416 18141 : libmesh_assert_less (ns, neighbor->n_neighbors());
417 :
418 : // Get all the active children (& grandchildren, etc...)
419 : // of the neighbor.
420 98801 : neighbor->active_family_tree (neighbors_offspring);
421 :
422 : // Get all the neighbor's children that
423 : // live on that side and are thus connected
424 : // to us
425 584297 : for (const auto & child : neighbors_offspring)
426 : {
427 : // Skip neighbor offspring which are not in the range of elements being partitioned.
428 485496 : if (!global_index_map.count(child->id()))
429 0 : continue;
430 :
431 : // This does not assume a level-1 mesh.
432 : // Note that since children have sides numbered
433 : // coincident with the parent then this is a sufficient test.
434 579031 : if (child->neighbor_ptr(ns) == elem)
435 : {
436 40870 : libmesh_assert (child->active());
437 :
438 211079 : csr_graph(elem_global_index, connection++) = global_index_map[child->id()];
439 : }
440 : }
441 : }
442 :
443 : #endif /* ifdef LIBMESH_ENABLE_AMR */
444 :
445 : }
446 : }
447 :
448 3913240 : if ((elem->dim() < LIBMESH_DIM) &&
449 1371100 : elem->interior_parent())
450 : {
451 : // get all relevant interior elements
452 440 : std::set<const Elem *> neighbor_set;
453 877 : elem->find_interior_neighbors(neighbor_set);
454 :
455 2064 : for (const Elem * neighbor : neighbor_set)
456 : {
457 : // Not all interior neighbors are necessarily in
458 : // the same Mesh (hence not in the global_index_map).
459 : // This will be the case when partitioning a
460 : // BoundaryMesh, whose elements all have
461 : // interior_parents() that belong to some other
462 : // Mesh.
463 1187 : const Elem * queried_elem = mesh.query_elem_ptr(neighbor->id());
464 :
465 : // Compare the neighbor and the queried_elem
466 : // pointers, make sure they are the same.
467 1187 : if (queried_elem && queried_elem == neighbor)
468 1689 : csr_graph(elem_global_index, connection++) =
469 1187 : libmesh_map_find(global_index_map, neighbor->id());
470 : }
471 : }
472 :
473 : // Check for any boundary neighbors
474 2543327 : for (const auto & pr : as_range(interior_to_boundary_map.equal_range(elem)))
475 : {
476 1187 : const Elem * neighbor = pr.second;
477 1438 : csr_graph(elem_global_index, connection++) =
478 1187 : global_index_map[neighbor->id()];
479 : }
480 25961 : }
481 :
482 : // We create a non-empty vals for a disconnected graph, to
483 : // work around a segfault from METIS.
484 4000 : libmesh_assert_equal_to (csr_graph.vals.size(),
485 : std::max(graph_size, std::size_t(1)));
486 : } // done building the graph
487 :
488 33961 : Metis::idx_t ncon = 1;
489 :
490 : // Select which type of partitioning to create
491 :
492 : // Use recursive if the number of partitions is less than or equal to 8
493 33961 : if (n_pieces <= 8)
494 28369 : Metis::METIS_PartGraphRecursive(&n,
495 : &ncon,
496 : csr_graph.offsets.data(),
497 : csr_graph.vals.data(),
498 : vwgt.data(),
499 : nullptr,
500 : nullptr,
501 : &nparts,
502 : nullptr,
503 : nullptr,
504 : nullptr,
505 : &edgecut,
506 : part.data());
507 :
508 : // Otherwise use kway
509 : else
510 5592 : Metis::METIS_PartGraphKway(&n,
511 : &ncon,
512 : csr_graph.offsets.data(),
513 : csr_graph.vals.data(),
514 : vwgt.data(),
515 : nullptr,
516 : nullptr,
517 : &nparts,
518 : nullptr,
519 : nullptr,
520 : nullptr,
521 : &edgecut,
522 : part.data());
523 :
524 : } // end processor 0 part
525 :
526 : // Broadcast the resulting partition
527 233031 : mesh.comm().broadcast(part);
528 :
529 : // Assign the returned processor ids. The part array contains
530 : // the processor id for each active element, but in terms of
531 : // the contiguous indexing we defined above
532 13322829 : for (auto & elem : as_range(beg, end))
533 : {
534 1046528 : libmesh_assert (global_index_map.count(elem->id()));
535 :
536 : const dof_id_type elem_global_index =
537 12864767 : global_index_map[elem->id()];
538 :
539 1046528 : libmesh_assert_less (elem_global_index, part.size());
540 : const processor_id_type elem_procid =
541 12864767 : static_cast<processor_id_type>(part[elem_global_index]);
542 :
543 12864767 : elem->processor_id() = elem_procid;
544 217031 : }
545 : #endif
546 : }
547 :
548 :
549 :
550 233098 : void MetisPartitioner::_do_partition (MeshBase & mesh,
551 : const unsigned int n_pieces)
552 : {
553 233098 : this->partition_range(mesh,
554 466196 : mesh.active_elements_begin(),
555 241100 : mesh.active_elements_end(),
556 16004 : n_pieces);
557 233098 : }
558 :
559 : } // namespace libMesh
|
