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