Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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/parmetis_partitioner.h"
22 :
23 : // libMesh includes
24 : #include "libmesh/libmesh_config.h"
25 : #include "libmesh/elem.h"
26 : #include "libmesh/elem_range.h"
27 : #include "libmesh/enum_partitioner_type.h"
28 : #include "libmesh/libmesh_logging.h"
29 : #include "libmesh/mesh_base.h"
30 : #include "libmesh/mesh_serializer.h"
31 : #include "libmesh/mesh_tools.h"
32 : #include "libmesh/mesh_communication.h"
33 : #include "libmesh/metis_partitioner.h"
34 : #include "libmesh/parallel_only.h"
35 : #include "libmesh/parmetis_helper.h"
36 : #include "libmesh/utility.h"
37 :
38 : // TIMPI includes
39 : #include "timpi/communicator.h" // also includes mpi.h
40 : #include "timpi/parallel_implementation.h" // for min()
41 :
42 : // Include the ParMETIS header file.
43 : #ifdef LIBMESH_HAVE_PARMETIS
44 :
45 : // Before we include a header wrapped in a namespace, we'd better make
46 : // sure none of its dependencies end up in that namespace
47 : #include <mpi.h>
48 :
49 : namespace Parmetis {
50 : extern "C" {
51 : # include "libmesh/ignore_warnings.h"
52 : # include "parmetis.h"
53 : # include "libmesh/restore_warnings.h"
54 : }
55 : }
56 :
57 : #endif
58 :
59 :
60 : // C++ includes
61 : #include <unordered_map>
62 : #include <unordered_set>
63 :
64 :
65 : namespace libMesh
66 : {
67 :
68 : // Minimum elements on each processor required for us to choose
69 : // Parmetis over Metis.
70 : #ifdef LIBMESH_HAVE_PARMETIS
71 : const unsigned int MIN_ELEM_PER_PROC = 4;
72 : #endif
73 :
74 : // ------------------------------------------------------------
75 : // ParmetisPartitioner implementation
76 278660 : ParmetisPartitioner::ParmetisPartitioner()
77 : #ifdef LIBMESH_HAVE_PARMETIS
78 278932 : : _pmetis(std::make_unique<ParmetisHelper>())
79 : #endif
80 278660 : {}
81 :
82 :
83 :
84 19326 : ParmetisPartitioner::ParmetisPartitioner (const ParmetisPartitioner & other)
85 : : Partitioner(other)
86 : #ifdef LIBMESH_HAVE_PARMETIS
87 19326 : , _pmetis(std::make_unique<ParmetisHelper>(*(other._pmetis)))
88 : #endif
89 : {
90 19326 : }
91 :
92 :
93 :
94 594070 : ParmetisPartitioner::~ParmetisPartitioner() = default;
95 :
96 :
97 :
98 284 : PartitionerType ParmetisPartitioner::type() const
99 : {
100 284 : return PARMETIS_PARTITIONER;
101 : }
102 :
103 :
104 :
105 260030 : void ParmetisPartitioner::_do_partition (MeshBase & mesh,
106 : const unsigned int n_sbdmns)
107 : {
108 260030 : this->_do_repartition (mesh, n_sbdmns);
109 260030 : }
110 :
111 :
112 :
113 260030 : void ParmetisPartitioner::_do_repartition (MeshBase & mesh,
114 : const unsigned int n_sbdmns)
115 : {
116 : // This function must be run on all processors at once
117 380 : libmesh_parallel_only(mesh.comm());
118 :
119 : // Check for easy returns
120 260030 : if (!mesh.n_elem())
121 211036 : return;
122 :
123 259888 : if (n_sbdmns == 1)
124 : {
125 0 : this->single_partition(mesh);
126 0 : return;
127 : }
128 :
129 376 : libmesh_assert_greater (n_sbdmns, 0);
130 :
131 : // What to do if the Parmetis library IS NOT present
132 : #ifndef LIBMESH_HAVE_PARMETIS
133 :
134 : libmesh_do_once(
135 : libMesh::out << "ERROR: The library has been built without" << std::endl
136 : << "Parmetis support. Using a Metis" << std::endl
137 : << "partitioner instead!" << std::endl;);
138 :
139 : MetisPartitioner mp;
140 :
141 : // Metis and other fallbacks only work in serial, and need to get
142 : // handed element ranges from an already-serialized mesh.
143 : mesh.allgather();
144 :
145 : // Don't just call partition() here; that would end up calling
146 : // post-element-partitioning work redundantly (and at the moment
147 : // incorrectly)
148 : mp.partition_range (mesh, mesh.active_elements_begin(),
149 : mesh.active_elements_end(), n_sbdmns);
150 :
151 : // What to do if the Parmetis library IS present
152 : #else
153 :
154 : // Revert to METIS on one processor.
155 260264 : if (mesh.n_processors() == 1)
156 : {
157 : // Make sure the mesh knows it's serial
158 4 : mesh.allgather();
159 :
160 0 : MetisPartitioner mp;
161 : // Don't just call partition() here; that would end up calling
162 : // post-element-partitioning work redundantly (and at the moment
163 : // incorrectly)
164 8 : mp.partition_range (mesh, mesh.active_elements_begin(),
165 8 : mesh.active_elements_end(), n_sbdmns);
166 0 : return;
167 : }
168 :
169 376 : LOG_SCOPE("repartition()", "ParmetisPartitioner");
170 :
171 : // Initialize the data structures required by ParMETIS
172 259884 : this->initialize (mesh, n_sbdmns);
173 :
174 : // build the graph corresponding to the mesh
175 259884 : this->build_graph (mesh);
176 :
177 : // Make sure all processors have enough active local elements and
178 : // enough connectivity among them.
179 : // Parmetis tends to die when it's given only a couple elements
180 : // per partition or when it can't reach elements from each other.
181 : {
182 376 : bool ready_for_parmetis = true;
183 3050876 : for (const auto & nelem : _n_active_elem_on_proc)
184 2790992 : if (nelem < MIN_ELEM_PER_PROC)
185 154 : ready_for_parmetis = false;
186 :
187 259884 : std::size_t my_adjacency = _pmetis->adjncy.size();
188 259884 : mesh.comm().min(my_adjacency);
189 259884 : if (!my_adjacency)
190 18 : ready_for_parmetis = false;
191 :
192 : // Parmetis will not work unless each processor has some
193 : // elements. Specifically, it will abort when passed a nullptr
194 : // partition or adjacency array on *any* of the processors.
195 107391 : if (!ready_for_parmetis)
196 : {
197 : // FIXME: revert to METIS, although this requires a serial mesh
198 211070 : MeshSerializer serialize(mesh);
199 90 : MetisPartitioner mp;
200 210890 : mp.partition (mesh, n_sbdmns);
201 90 : return;
202 210710 : }
203 : }
204 :
205 :
206 : // Partition the graph
207 49852 : std::vector<Parmetis::idx_t> vsize(_pmetis->vwgt.size(), 1);
208 48994 : Parmetis::real_t itr = 1000000.0;
209 48994 : MPI_Comm mpi_comm = mesh.comm().get();
210 :
211 : // Call the ParMETIS adaptive repartitioning method. This respects the
212 : // original partitioning when computing the new partitioning so as to
213 : // minimize the required data redistribution.
214 441518 : Parmetis::ParMETIS_V3_AdaptiveRepart(_pmetis->vtxdist.empty() ? nullptr : _pmetis->vtxdist.data(),
215 572 : _pmetis->xadj.empty() ? nullptr : _pmetis->xadj.data(),
216 572 : _pmetis->adjncy.empty() ? nullptr : _pmetis->adjncy.data(),
217 572 : _pmetis->vwgt.empty() ? nullptr : _pmetis->vwgt.data(),
218 572 : vsize.empty() ? nullptr : vsize.data(),
219 : nullptr,
220 286 : &_pmetis->wgtflag,
221 286 : &_pmetis->numflag,
222 286 : &_pmetis->ncon,
223 286 : &_pmetis->nparts,
224 572 : _pmetis->tpwgts.empty() ? nullptr : _pmetis->tpwgts.data(),
225 572 : _pmetis->ubvec.empty() ? nullptr : _pmetis->ubvec.data(),
226 : &itr,
227 286 : _pmetis->options.data(),
228 286 : &_pmetis->edgecut,
229 572 : _pmetis->part.empty() ? nullptr : reinterpret_cast<Parmetis::idx_t *>(_pmetis->part.data()),
230 : &mpi_comm);
231 :
232 : // Assign the returned processor ids
233 49280 : this->assign_partitioning (mesh, _pmetis->part);
234 :
235 : #endif // #ifndef LIBMESH_HAVE_PARMETIS ... else ...
236 :
237 : }
238 :
239 :
240 :
241 : // Only need to compile these methods if ParMETIS is present
242 : #ifdef LIBMESH_HAVE_PARMETIS
243 :
244 259884 : void ParmetisPartitioner::initialize (const MeshBase & mesh,
245 : const unsigned int n_sbdmns)
246 : {
247 752 : LOG_SCOPE("initialize()", "ParmetisPartitioner");
248 :
249 376 : const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
250 : // Set parameters.
251 259884 : _pmetis->wgtflag = 2; // weights on vertices only
252 259884 : _pmetis->ncon = 1; // one weight per vertex
253 259884 : _pmetis->numflag = 0; // C-style 0-based numbering
254 259884 : _pmetis->nparts = static_cast<Parmetis::idx_t>(n_sbdmns); // number of subdomains to create
255 259884 : _pmetis->edgecut = 0; // the numbers of edges cut by the
256 : // partition
257 :
258 : // Initialize data structures for ParMETIS
259 260260 : _pmetis->vtxdist.assign (mesh.n_processors()+1, 0);
260 259884 : _pmetis->tpwgts.assign (_pmetis->nparts, 1./_pmetis->nparts);
261 259884 : _pmetis->ubvec.assign (_pmetis->ncon, 1.05);
262 259884 : _pmetis->part.assign (n_active_local_elem, 0);
263 259884 : _pmetis->options.resize (5);
264 259884 : _pmetis->vwgt.resize (n_active_local_elem);
265 :
266 : // Set the options
267 259884 : _pmetis->options[0] = 1; // don't use default options
268 259884 : _pmetis->options[1] = 0; // default (level of timing)
269 259884 : _pmetis->options[2] = 15; // random seed (default)
270 259884 : _pmetis->options[3] = 2; // processor distribution and subdomain distribution are decoupled
271 :
272 : // ParMetis expects the elements to be numbered in contiguous blocks
273 : // by processor, i.e. [0, ne0), [ne0, ne0+ne1), ...
274 : // Since we only partition active elements we should have no expectation
275 : // that we currently have such a distribution. So we need to create it.
276 : // Also, at the same time we are going to map all the active elements into a globally
277 : // unique range [0,n_active_elem) which is *independent* of the current partitioning.
278 : // This can be fed to ParMetis as the initial partitioning of the subdomains (decoupled
279 : // from the partitioning of the objects themselves). This allows us to get the same
280 : // resultant partitioning independent of the input partitioning.
281 : libMesh::BoundingBox bbox =
282 259884 : MeshTools::create_bounding_box(mesh);
283 :
284 259884 : _find_global_index_by_pid_map(mesh);
285 :
286 :
287 : // count the total number of active elements in the mesh. Note we cannot
288 : // use mesh.n_active_elem() in general since this only returns the number
289 : // of active elements which are stored on the calling processor.
290 : // We should not use n_active_elem for any allocation because that will
291 : // be inherently unscalable, but it can be useful for libmesh_assertions.
292 376 : dof_id_type n_active_elem=0;
293 :
294 : // Set up the vtxdist array. This will be the same on each processor.
295 : // ***** Consult the Parmetis documentation. *****
296 376 : libmesh_assert_equal_to (_pmetis->vtxdist.size(),
297 : cast_int<std::size_t>(mesh.n_processors()+1));
298 376 : libmesh_assert_equal_to (_pmetis->vtxdist[0], 0);
299 :
300 3050876 : for (auto pid : make_range(mesh.n_processors()))
301 : {
302 2792496 : _pmetis->vtxdist[pid+1] = _pmetis->vtxdist[pid] + _n_active_elem_on_proc[pid];
303 2790992 : n_active_elem += _n_active_elem_on_proc[pid];
304 : }
305 376 : libmesh_assert_equal_to (_pmetis->vtxdist.back(), static_cast<Parmetis::idx_t>(n_active_elem));
306 :
307 :
308 : // Maps active element ids into a contiguous range independent of partitioning.
309 : // (only needs local scope)
310 752 : std::unordered_map<dof_id_type, dof_id_type> global_index_map;
311 :
312 : {
313 752 : std::vector<dof_id_type> global_index;
314 :
315 : // create the unique mapping for all active elements independent of partitioning
316 : {
317 260260 : MeshBase::const_element_iterator it = mesh.active_elements_begin();
318 260260 : const MeshBase::const_element_iterator end = mesh.active_elements_end();
319 :
320 : // Calling this on all processors a unique range in [0,n_active_elem) is constructed.
321 : // Only the indices for the elements we pass in are returned in the array.
322 259884 : MeshCommunication().find_global_indices (mesh.comm(),
323 : bbox, it, end,
324 : global_index);
325 :
326 16816446 : for (dof_id_type cnt=0; it != end; ++it)
327 : {
328 16297430 : const Elem * elem = *it;
329 : // vectormap::count forces a sort, which is too expensive
330 : // in a loop
331 : // libmesh_assert (!global_index_map.count(elem->id()));
332 19674 : libmesh_assert_less (cnt, global_index.size());
333 19674 : libmesh_assert_less (global_index[cnt], n_active_elem);
334 :
335 16317104 : global_index_map.emplace(elem->id(), global_index[cnt++]);
336 : }
337 : }
338 : // really, shouldn't be close!
339 376 : libmesh_assert_less_equal (global_index_map.size(), n_active_elem);
340 376 : libmesh_assert_less_equal (_global_index_by_pid_map.size(), n_active_elem);
341 :
342 : // At this point the two maps should be the same size. If they are not
343 : // then the number of active elements is not the same as the sum over all
344 : // processors of the number of active elements per processor, which means
345 : // there must be some unpartitioned objects out there.
346 259884 : libmesh_error_msg_if(global_index_map.size() != _global_index_by_pid_map.size(),
347 : "ERROR: ParmetisPartitioner cannot handle unpartitioned objects!");
348 : }
349 :
350 : // Finally, we need to initialize the vertex (partition) weights and the initial subdomain
351 : // mapping. The subdomain mapping will be independent of the processor mapping, and is
352 : // defined by a simple mapping of the global indices we just found.
353 : {
354 260260 : std::vector<dof_id_type> subdomain_bounds(mesh.n_processors());
355 :
356 260636 : const dof_id_type first_local_elem = _pmetis->vtxdist[mesh.processor_id()];
357 :
358 3050876 : for (auto pid : make_range(mesh.n_processors()))
359 : {
360 752 : dof_id_type tgt_subdomain_size = 0;
361 :
362 : // watch out for the case that n_subdomains < n_processors
363 2790992 : if (pid < static_cast<unsigned int>(_pmetis->nparts))
364 : {
365 1797664 : tgt_subdomain_size = n_active_elem/std::min
366 1797664 : (cast_int<Parmetis::idx_t>(mesh.n_processors()), _pmetis->nparts);
367 :
368 1796912 : if (pid < n_active_elem%_pmetis->nparts)
369 620353 : tgt_subdomain_size++;
370 : }
371 2790992 : if (pid == 0)
372 259884 : subdomain_bounds[0] = tgt_subdomain_size;
373 : else
374 2531860 : subdomain_bounds[pid] = subdomain_bounds[pid-1] + tgt_subdomain_size;
375 : }
376 :
377 376 : libmesh_assert_equal_to (subdomain_bounds.back(), n_active_elem);
378 :
379 : Threads::parallel_for
380 259884 : (mesh.active_local_element_stored_range(),
381 777396 : [first_local_elem, n_active_elem, n_active_local_elem, this,
382 166603 : &global_index_map, &subdomain_bounds](const ConstElemRange & range)
383 : {
384 2940905 : for (const Elem * elem : range)
385 : {
386 : const dof_id_type global_index_by_pid =
387 2681021 : libmesh_map_find(_global_index_by_pid_map, elem->id());
388 11774 : libmesh_assert_less (global_index_by_pid, n_active_elem);
389 :
390 2681021 : const dof_id_type local_index =
391 23548 : global_index_by_pid - first_local_elem;
392 :
393 11774 : libmesh_assert_less (local_index, n_active_local_elem);
394 11774 : libmesh_assert_less (local_index, _pmetis->vwgt.size());
395 :
396 11774 : libmesh_ignore(n_active_elem); // unused outside dbg/devel
397 11774 : libmesh_ignore(n_active_local_elem); // unused outside dbg/devel
398 :
399 : // Spline nodes are a special case (storing all the
400 : // unconstrained DoFs in an IGA simulation), but in general
401 : // we'll try to distribute work by expecting it to be roughly
402 : // proportional to DoFs, which are roughly proportional to
403 : // nodes.
404 2681021 : if (elem->type() == NODEELEM &&
405 0 : elem->mapping_type() == RATIONAL_BERNSTEIN_MAP)
406 12789 : _pmetis->vwgt[local_index] = 50;
407 : else
408 2668232 : _pmetis->vwgt[local_index] = elem->n_nodes();
409 :
410 : // find the subdomain this element belongs in
411 : const dof_id_type global_index =
412 2681021 : libmesh_map_find(global_index_map, elem->id());
413 :
414 11774 : libmesh_assert_less (global_index, subdomain_bounds.back());
415 :
416 : const unsigned int subdomain_id =
417 : cast_int<unsigned int>
418 2669247 : (std::distance(subdomain_bounds.begin(),
419 : std::lower_bound(subdomain_bounds.begin(),
420 : subdomain_bounds.end(),
421 11774 : global_index)));
422 11774 : libmesh_assert_less (subdomain_id, _pmetis->nparts);
423 11774 : libmesh_assert_less (local_index, _pmetis->part.size());
424 :
425 2692795 : _pmetis->part[local_index] = subdomain_id;
426 : }
427 259884 : });
428 : }
429 259884 : }
430 :
431 :
432 :
433 259884 : void ParmetisPartitioner::build_graph (const MeshBase & mesh)
434 : {
435 752 : LOG_SCOPE("build_graph()", "ParmetisPartitioner");
436 :
437 : // build the graph in distributed CSR format. Note that
438 : // the edges in the graph will correspond to
439 : // face neighbors
440 376 : const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
441 :
442 259884 : Partitioner::build_graph(mesh);
443 :
444 376 : dof_id_type graph_size_upper_bound = 0;
445 :
446 2940905 : for (auto & row: _dual_graph)
447 2692795 : graph_size_upper_bound += cast_int<dof_id_type>(row.size());
448 :
449 : // Reserve space in the adjacency array
450 259884 : _pmetis->xadj.clear();
451 259884 : _pmetis->xadj.reserve (n_active_local_elem + 1);
452 259884 : _pmetis->adjncy.clear();
453 259884 : _pmetis->adjncy.reserve (graph_size_upper_bound);
454 :
455 : const dof_id_type first_local_elem =
456 260636 : _pmetis->vtxdist[mesh.processor_id()];
457 752 : std::unordered_set<dof_id_type> unique_neighbors;
458 :
459 : // ParMETIS requires a simple graph. Preserve first-occurrence order
460 : // while filtering so valid graph rows remain otherwise unchanged.
461 2940905 : for (auto i : index_range(_dual_graph))
462 : {
463 2692795 : _pmetis->xadj.push_back(cast_int<int>(_pmetis->adjncy.size()));
464 :
465 : const dof_id_type global_index =
466 2681021 : first_local_elem + cast_int<dof_id_type>(i);
467 23548 : const auto & graph_row = _dual_graph[i];
468 11774 : unique_neighbors.clear();
469 :
470 13653451 : for (const auto neighbor : graph_row)
471 11023308 : if (neighbor != global_index &&
472 10972052 : unique_neighbors.insert(neighbor).second)
473 10970988 : _pmetis->adjncy.push_back(neighbor);
474 : }
475 :
476 : // The end of the adjacency array for the last elem
477 260260 : _pmetis->xadj.push_back(cast_int<int>(_pmetis->adjncy.size()));
478 :
479 376 : libmesh_assert_equal_to (_pmetis->xadj.size(), n_active_local_elem+1);
480 376 : libmesh_assert_less_equal (_pmetis->adjncy.size(),
481 : graph_size_upper_bound);
482 259884 : }
483 :
484 : #endif // #ifdef LIBMESH_HAVE_PARMETIS
485 :
486 : } // namespace libMesh
|
