metis_partitioner.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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/libmesh_config.h"
22 #include "libmesh/mesh_base.h"
25 #include "libmesh/elem.h"
27 #include "libmesh/error_vector.h"
28 #include "libmesh/vectormap.h"
30 
31 #ifdef LIBMESH_HAVE_METIS
32 // MIPSPro 7.4.2 gets confused about these nested namespaces
33 # ifdef __sgi
34 # include <cstdarg>
35 # endif
36 namespace Metis {
37 extern "C" {
38 # include "libmesh/ignore_warnings.h"
39 # include "metis.h"
41 }
42 }
43 #else
44 # include "libmesh/sfc_partitioner.h"
45 #endif
46 
47 
48 // C++ includes
49 #include <unordered_map>
50 
51 
52 namespace libMesh
53 {
54 
55 
59  unsigned int n_pieces)
60 {
61  // Check for easy returns
62  if (beg == end)
63  return;
64 
65  if (n_pieces == 1)
66  {
67  this->single_partition_range (beg, end);
68  return;
69  }
70 
71  libmesh_assert_greater (n_pieces, 0);
72 
73  // We don't yet support distributed meshes with this Partitioner
74  if (!mesh.is_serial())
75  {
76  libMesh::out << "WARNING: Forced to gather a distributed mesh for METIS" << std::endl;
77  mesh.allgather();
78  }
79 
80  // What to do if the Metis library IS NOT present
81 #ifndef LIBMESH_HAVE_METIS
82 
83  libmesh_do_once(
84  libMesh::out << "ERROR: The library has been built without" << std::endl
85  << "Metis support. Using a space-filling curve" << std::endl
86  << "partitioner instead!" << std::endl;);
87 
88  SFCPartitioner sfcp;
89  sfcp.partition_range (mesh, beg, end, n_pieces);
90 
91  // What to do if the Metis library IS present
92 #else
93 
94  LOG_SCOPE("partition_range()", "MetisPartitioner");
95 
96  const dof_id_type n_range_elem = cast_int<dof_id_type>(std::distance(beg, end));
97 
98  // Metis will only consider the elements in the range.
99  // We need to map the range element ids into a
100  // contiguous range. Further, we want the unique range indexing to be
101  // independent of the element ordering, otherwise a circular dependency
102  // can result in which the partitioning depends on the ordering which
103  // depends on the partitioning...
104  vectormap<dof_id_type, dof_id_type> global_index_map;
105  global_index_map.reserve (n_range_elem);
106 
107  {
108  std::vector<dof_id_type> global_index;
109 
112  beg, end, global_index);
113 
114  libmesh_assert_equal_to (global_index.size(), n_range_elem);
115 
116  std::size_t cnt=0;
117  for (const auto & elem : as_range(beg, end))
118  global_index_map.insert (std::make_pair(elem->id(), global_index[cnt++]));
119 
120  libmesh_assert_equal_to (global_index_map.size(), n_range_elem);
121  }
122 
123  // If we have boundary elements in this mesh, we want to account for
124  // the connectivity between them and interior elements. We can find
125  // interior elements from boundary elements, but we need to build up
126  // a lookup map to do the reverse.
127  typedef std::unordered_multimap<const Elem *, const Elem *> map_type;
128  map_type interior_to_boundary_map;
129 
130  for (const auto & elem : as_range(beg, end))
131  {
132  // If we don't have an interior_parent then there's nothing
133  // to look us up.
134  if ((elem->dim() >= LIBMESH_DIM) ||
135  !elem->interior_parent())
136  continue;
137 
138  // get all relevant interior elements
139  std::set<const Elem *> neighbor_set;
140  elem->find_interior_neighbors(neighbor_set);
141 
142  std::set<const Elem *>::iterator n_it = neighbor_set.begin();
143  for (; n_it != neighbor_set.end(); ++n_it)
144  {
145  // FIXME - non-const versions of the std::set<const Elem
146  // *> returning methods would be nice
147  Elem * neighbor = const_cast<Elem *>(*n_it);
148 
149 #if defined(LIBMESH_HAVE_UNORDERED_MULTIMAP) || \
150  defined(LIBMESH_HAVE_TR1_UNORDERED_MULTIMAP) || \
151  defined(LIBMESH_HAVE_HASH_MULTIMAP) || \
152  defined(LIBMESH_HAVE_EXT_HASH_MULTIMAP)
153  interior_to_boundary_map.insert(std::make_pair(neighbor, elem));
154 #else
155  interior_to_boundary_map.insert(interior_to_boundary_map.begin(),
156  std::make_pair(neighbor, elem));
157 #endif
158  }
159  }
160 
161  // Data structure that Metis will fill up on processor 0 and broadcast.
162  std::vector<Metis::idx_t> part(n_range_elem);
163 
164  // Invoke METIS, but only on processor 0.
165  // Then broadcast the resulting decomposition
166  if (mesh.processor_id() == 0)
167  {
168  // Data structures and parameters needed only on processor 0 by Metis.
169  // std::vector<Metis::idx_t> options(5);
170  std::vector<Metis::idx_t> vwgt(n_range_elem);
171 
172  Metis::idx_t
173  n = static_cast<Metis::idx_t>(n_range_elem), // number of "nodes" (elements) in the graph
174  // wgtflag = 2, // weights on vertices only, none on edges
175  // numflag = 0, // C-style 0-based numbering
176  nparts = static_cast<Metis::idx_t>(n_pieces), // number of subdomains to create
177  edgecut = 0; // the numbers of edges cut by the resulting partition
178 
179  // Set the options
180  // options[0] = 0; // use default options
181 
182  // build the graph
184 
185  csr_graph.offsets.resize(n_range_elem + 1, 0);
186 
187  // Local scope for these
188  {
189  // build the graph in CSR format. Note that
190  // the edges in the graph will correspond to
191  // face neighbors
192 
193 #ifdef LIBMESH_ENABLE_AMR
194  std::vector<const Elem *> neighbors_offspring;
195 #endif
196 
197 #ifndef NDEBUG
198  std::size_t graph_size=0;
199 #endif
200 
201  // (1) first pass - get the row sizes for each element by counting the number
202  // of face neighbors. Also populate the vwght array if necessary
203  for (const auto & elem : as_range(beg, end))
204  {
205  const dof_id_type elem_global_index =
206  global_index_map[elem->id()];
207 
208  libmesh_assert_less (elem_global_index, vwgt.size());
209 
210  // maybe there is a better weight?
211  // The weight is used to define what a balanced graph is
212  if (!_weights)
213  vwgt[elem_global_index] = elem->n_nodes();
214  else
215  vwgt[elem_global_index] = static_cast<Metis::idx_t>((*_weights)[elem->id()]);
216 
217  unsigned int num_neighbors = 0;
218 
219  // Loop over the element's neighbors. An element
220  // adjacency corresponds to a face neighbor
221  for (auto neighbor : elem->neighbor_ptr_range())
222  {
223  if (neighbor != nullptr)
224  {
225  // If the neighbor is active, but is not in the
226  // range of elements being partitioned, treat it
227  // as a nullptr neighbor.
228  if (neighbor->active() && !global_index_map.count(neighbor->id()))
229  continue;
230 
231  // If the neighbor is active treat it
232  // as a connection
233  if (neighbor->active())
234  num_neighbors++;
235 
236 #ifdef LIBMESH_ENABLE_AMR
237 
238  // Otherwise we need to find all of the
239  // neighbor's children that are connected to
240  // us and add them
241  else
242  {
243  // The side of the neighbor to which
244  // we are connected
245  const unsigned int ns =
246  neighbor->which_neighbor_am_i (elem);
247  libmesh_assert_less (ns, neighbor->n_neighbors());
248 
249  // Get all the active children (& grandchildren, etc...)
250  // of the neighbor.
251 
252  // FIXME - this is the wrong thing, since we
253  // should be getting the active family tree on
254  // our side only. But adding too many graph
255  // links may cause hanging nodes to tend to be
256  // on partition interiors, which would reduce
257  // communication overhead for constraint
258  // equations, so we'll leave it.
259  neighbor->active_family_tree (neighbors_offspring);
260 
261  // Get all the neighbor's children that
262  // live on that side and are thus connected
263  // to us
264  for (std::size_t nc=0; nc<neighbors_offspring.size(); nc++)
265  {
266  const Elem * child =
267  neighbors_offspring[nc];
268 
269  // Skip neighbor offspring which are not in the range of elements being partitioned.
270  if (!global_index_map.count(child->id()))
271  continue;
272 
273  // This does not assume a level-1 mesh.
274  // Note that since children have sides numbered
275  // coincident with the parent then this is a sufficient test.
276  if (child->neighbor_ptr(ns) == elem)
277  {
278  libmesh_assert (child->active());
279  num_neighbors++;
280  }
281  }
282  }
283 
284 #endif /* ifdef LIBMESH_ENABLE_AMR */
285 
286  }
287  }
288 
289  // Check for any interior neighbors
290  if ((elem->dim() < LIBMESH_DIM) && elem->interior_parent())
291  {
292  // get all relevant interior elements
293  std::set<const Elem *> neighbor_set;
294  elem->find_interior_neighbors(neighbor_set);
295 
296  num_neighbors += cast_int<unsigned int>(neighbor_set.size());
297  }
298 
299  // Check for any boundary neighbors
300  typedef map_type::iterator map_it_type;
301  std::pair<map_it_type, map_it_type>
302  bounds = interior_to_boundary_map.equal_range(elem);
303  num_neighbors += cast_int<unsigned int>
304  (std::distance(bounds.first, bounds.second));
305 
306  csr_graph.prep_n_nonzeros(elem_global_index, num_neighbors);
307 #ifndef NDEBUG
308  graph_size += num_neighbors;
309 #endif
310  }
311 
312  csr_graph.prepare_for_use();
313 
314  // (2) second pass - fill the compressed adjacency array
315  for (const auto & elem : as_range(beg, end))
316  {
317  const dof_id_type elem_global_index =
318  global_index_map[elem->id()];
319 
320  unsigned int connection=0;
321 
322  // Loop over the element's neighbors. An element
323  // adjacency corresponds to a face neighbor
324  for (auto neighbor : elem->neighbor_ptr_range())
325  {
326  if (neighbor != nullptr)
327  {
328  // If the neighbor is active, but is not in the
329  // range of elements being partitioned, treat it
330  // as a nullptr neighbor.
331  if (neighbor->active() && !global_index_map.count(neighbor->id()))
332  continue;
333 
334  // If the neighbor is active treat it
335  // as a connection
336  if (neighbor->active())
337  csr_graph(elem_global_index, connection++) = global_index_map[neighbor->id()];
338 
339 #ifdef LIBMESH_ENABLE_AMR
340 
341  // Otherwise we need to find all of the
342  // neighbor's children that are connected to
343  // us and add them
344  else
345  {
346  // The side of the neighbor to which
347  // we are connected
348  const unsigned int ns =
349  neighbor->which_neighbor_am_i (elem);
350  libmesh_assert_less (ns, neighbor->n_neighbors());
351 
352  // Get all the active children (& grandchildren, etc...)
353  // of the neighbor.
354  neighbor->active_family_tree (neighbors_offspring);
355 
356  // Get all the neighbor's children that
357  // live on that side and are thus connected
358  // to us
359  for (std::size_t nc=0; nc<neighbors_offspring.size(); nc++)
360  {
361  const Elem * child =
362  neighbors_offspring[nc];
363 
364  // Skip neighbor offspring which are not in the range of elements being partitioned.
365  if (!global_index_map.count(child->id()))
366  continue;
367 
368  // This does not assume a level-1 mesh.
369  // Note that since children have sides numbered
370  // coincident with the parent then this is a sufficient test.
371  if (child->neighbor_ptr(ns) == elem)
372  {
373  libmesh_assert (child->active());
374 
375  csr_graph(elem_global_index, connection++) = global_index_map[child->id()];
376  }
377  }
378  }
379 
380 #endif /* ifdef LIBMESH_ENABLE_AMR */
381 
382  }
383  }
384 
385  if ((elem->dim() < LIBMESH_DIM) &&
386  elem->interior_parent())
387  {
388  // get all relevant interior elements
389  std::set<const Elem *> neighbor_set;
390  elem->find_interior_neighbors(neighbor_set);
391 
392  std::set<const Elem *>::iterator n_it = neighbor_set.begin();
393  for (; n_it != neighbor_set.end(); ++n_it)
394  {
395  const Elem * neighbor = *n_it;
396 
397  // Not all interior neighbors are necessarily in
398  // the same Mesh (hence not in the global_index_map).
399  // This will be the case when partitioning a
400  // BoundaryMesh, whose elements all have
401  // interior_parents() that belong to some other
402  // Mesh.
403  const Elem * queried_elem = mesh.query_elem_ptr(neighbor->id());
404 
405  // Compare the neighbor and the queried_elem
406  // pointers, make sure they are the same.
407  if (queried_elem && queried_elem == neighbor)
408  {
410  global_index_map.find(neighbor->id());
411 
412  // If the interior_neighbor is in the Mesh but
413  // not in the global_index_map, we have other issues.
414  if (global_index_map_it == global_index_map.end())
415  libmesh_error_msg("Interior neighbor with id " << neighbor->id() << " not found in global_index_map.");
416 
417  else
418  csr_graph(elem_global_index, connection++) = global_index_map_it->second;
419  }
420  }
421  }
422 
423  // Check for any boundary neighbors
424  for (const auto & pr : as_range(interior_to_boundary_map.equal_range(elem)))
425  {
426  const Elem * neighbor = pr.second;
427  csr_graph(elem_global_index, connection++) =
428  global_index_map[neighbor->id()];
429  }
430  }
431 
432  // We create a non-empty vals for a disconnected graph, to
433  // work around a segfault from METIS.
434  libmesh_assert_equal_to (csr_graph.vals.size(),
435  std::max(graph_size, std::size_t(1)));
436  } // done building the graph
437 
438  Metis::idx_t ncon = 1;
439 
440  // Select which type of partitioning to create
441 
442  // Use recursive if the number of partitions is less than or equal to 8
443  if (n_pieces <= 8)
444  Metis::METIS_PartGraphRecursive(&n,
445  &ncon,
446  csr_graph.offsets.data(),
447  csr_graph.vals.data(),
448  vwgt.data(),
449  nullptr,
450  nullptr,
451  &nparts,
452  nullptr,
453  nullptr,
454  nullptr,
455  &edgecut,
456  part.data());
457 
458  // Otherwise use kway
459  else
460  Metis::METIS_PartGraphKway(&n,
461  &ncon,
462  csr_graph.offsets.data(),
463  csr_graph.vals.data(),
464  vwgt.data(),
465  nullptr,
466  nullptr,
467  &nparts,
468  nullptr,
469  nullptr,
470  nullptr,
471  &edgecut,
472  part.data());
473 
474  } // end processor 0 part
475 
476  // Broadcast the resulting partition
477  mesh.comm().broadcast(part);
478 
479  // Assign the returned processor ids. The part array contains
480  // the processor id for each active element, but in terms of
481  // the contiguous indexing we defined above
482  for (auto & elem : as_range(beg, end))
483  {
484  libmesh_assert (global_index_map.count(elem->id()));
485 
486  const dof_id_type elem_global_index =
487  global_index_map[elem->id()];
488 
489  libmesh_assert_less (elem_global_index, part.size());
490  const processor_id_type elem_procid =
491  static_cast<processor_id_type>(part[elem_global_index]);
492 
493  elem->processor_id() = elem_procid;
494  }
495 #endif
496 }
497 
498 
499 
501  const unsigned int n_pieces)
502 {
503  this->partition_range(mesh,
504  mesh.active_elements_begin(),
505  mesh.active_elements_end(),
506  n_pieces);
507 }
508 
509 } // namespace libMesh
virtual void partition_range(MeshBase &mesh, MeshBase::element_iterator it, MeshBase::element_iterator end, const unsigned int n) override
void find_global_indices(const Parallel::Communicator &communicator, const libMesh::BoundingBox &, const ForwardIterator &, const ForwardIterator &, std::vector< dof_id_type > &) const
Compressed graph data structure used by MetisPartitioner.
libMesh::BoundingBox create_bounding_box(const MeshBase &mesh)
Definition: mesh_tools.C:386
void find_interior_neighbors(std::set< const Elem *> &neighbor_set) const
Definition: elem.C:725
virtual void partition_range(MeshBase &mesh, MeshBase::element_iterator it, MeshBase::element_iterator end, const unsigned int n) override
The base class for all geometric element types.
Definition: elem.h:100
MeshBase & mesh
uint8_t processor_id_type
Definition: id_types.h:99
IterBase * end
long double max(long double a, double b)
Base class for Mesh.
Definition: mesh_base.h:77
ErrorVector * _weights
Definition: partitioner.h:267
std::vector< IndexType > offsets
void single_partition_range(MeshBase::element_iterator it, MeshBase::element_iterator end)
Definition: partitioner.C:172
dof_id_type id() const
Definition: dof_object.h:655
SimpleRange< I > as_range(const std::pair< I, I > &p)
Definition: simple_range.h:57
std::vector< IndexType > vals
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2050
virtual void _do_partition(MeshBase &mesh, const unsigned int n) override
Partitioner based on different types of space filling curves.
vector_type::iterator iterator
Definition: vectormap.h:71
void prep_n_nonzeros(const libMesh::dof_id_type row, const libMesh::dof_id_type n_nonzeros_in)
OStreamProxy out(std::cout)
bool active() const
Definition: elem.h:2390
uint8_t dof_id_type
Definition: id_types.h:64