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 401269 : ParmetisPartitioner::ParmetisPartitioner()
77 : #ifdef LIBMESH_HAVE_PARMETIS
78 401555 : : _pmetis(std::make_unique<ParmetisHelper>())
79 : #endif
80 401269 : {}
81 :
82 :
83 :
84 19397 : ParmetisPartitioner::ParmetisPartitioner (const ParmetisPartitioner & other)
85 : : Partitioner(other)
86 : #ifdef LIBMESH_HAVE_PARMETIS
87 19397 : , _pmetis(std::make_unique<ParmetisHelper>(*(other._pmetis)))
88 : #endif
89 : {
90 19397 : }
91 :
92 :
93 :
94 839366 : ParmetisPartitioner::~ParmetisPartitioner() = default;
95 :
96 :
97 :
98 284 : PartitionerType ParmetisPartitioner::type() const
99 : {
100 284 : return PARMETIS_PARTITIONER;
101 : }
102 :
103 :
104 :
105 327419 : void ParmetisPartitioner::_do_partition (MeshBase & mesh,
106 : const unsigned int n_sbdmns)
107 : {
108 327419 : this->_do_repartition (mesh, n_sbdmns);
109 327419 : }
110 :
111 :
112 :
113 327419 : 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 394 : libmesh_parallel_only(mesh.comm());
118 :
119 : // Check for easy returns
120 327419 : if (!mesh.n_elem())
121 271845 : return;
122 :
123 327277 : if (n_sbdmns == 1)
124 : {
125 0 : this->single_partition(mesh);
126 0 : return;
127 : }
128 :
129 390 : 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 327667 : 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 684 : LOG_SCOPE("repartition()", "ParmetisPartitioner");
170 :
171 : // Initialize the data structures required by ParMETIS
172 327273 : this->initialize (mesh, n_sbdmns);
173 :
174 : // build the graph corresponding to the mesh
175 327273 : 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 390 : bool ready_for_parmetis = true;
183 3844200 : for (const auto & nelem : _n_active_elem_on_proc)
184 3516927 : if (nelem < MIN_ELEM_PER_PROC)
185 174 : ready_for_parmetis = false;
186 :
187 327273 : std::size_t my_adjacency = _pmetis->adjncy.size();
188 327273 : mesh.comm().min(my_adjacency);
189 327273 : if (!my_adjacency)
190 22 : 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 129504 : if (!ready_for_parmetis)
196 : {
197 : // FIXME: revert to METIS, although this requires a serial mesh
198 271891 : MeshSerializer serialize(mesh);
199 96 : MetisPartitioner mp;
200 271699 : mp.partition (mesh, n_sbdmns);
201 96 : return;
202 271507 : }
203 : }
204 :
205 :
206 : // Partition the graph
207 56456 : std::vector<Parmetis::idx_t> vsize(_pmetis->vwgt.size(), 1);
208 55574 : Parmetis::real_t itr = 1000000.0;
209 55574 : 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 500754 : Parmetis::ParMETIS_V3_AdaptiveRepart(_pmetis->vtxdist.empty() ? nullptr : _pmetis->vtxdist.data(),
215 588 : _pmetis->xadj.empty() ? nullptr : _pmetis->xadj.data(),
216 588 : _pmetis->adjncy.empty() ? nullptr : _pmetis->adjncy.data(),
217 588 : _pmetis->vwgt.empty() ? nullptr : _pmetis->vwgt.data(),
218 588 : vsize.empty() ? nullptr : vsize.data(),
219 : nullptr,
220 294 : &_pmetis->wgtflag,
221 294 : &_pmetis->numflag,
222 294 : &_pmetis->ncon,
223 294 : &_pmetis->nparts,
224 588 : _pmetis->tpwgts.empty() ? nullptr : _pmetis->tpwgts.data(),
225 588 : _pmetis->ubvec.empty() ? nullptr : _pmetis->ubvec.data(),
226 : &itr,
227 294 : _pmetis->options.data(),
228 294 : &_pmetis->edgecut,
229 588 : _pmetis->part.empty() ? nullptr : reinterpret_cast<Parmetis::idx_t *>(_pmetis->part.data()),
230 : &mpi_comm);
231 :
232 : // Assign the returned processor ids
233 55868 : 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 327273 : void ParmetisPartitioner::initialize (const MeshBase & mesh,
245 : const unsigned int n_sbdmns)
246 : {
247 1170 : LOG_SCOPE("initialize()", "ParmetisPartitioner");
248 :
249 390 : const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
250 : // Set parameters.
251 327273 : _pmetis->wgtflag = 2; // weights on vertices only
252 327273 : _pmetis->ncon = 1; // one weight per vertex
253 327273 : _pmetis->numflag = 0; // C-style 0-based numbering
254 327273 : _pmetis->nparts = static_cast<Parmetis::idx_t>(n_sbdmns); // number of subdomains to create
255 327273 : _pmetis->edgecut = 0; // the numbers of edges cut by the
256 : // partition
257 :
258 : // Initialize data structures for ParMETIS
259 327663 : _pmetis->vtxdist.assign (mesh.n_processors()+1, 0);
260 327273 : _pmetis->tpwgts.assign (_pmetis->nparts, 1./_pmetis->nparts);
261 327273 : _pmetis->ubvec.assign (_pmetis->ncon, 1.05);
262 327273 : _pmetis->part.assign (n_active_local_elem, 0);
263 327273 : _pmetis->options.resize (5);
264 327273 : _pmetis->vwgt.resize (n_active_local_elem);
265 :
266 : // Set the options
267 327273 : _pmetis->options[0] = 1; // don't use default options
268 327273 : _pmetis->options[1] = 0; // default (level of timing)
269 327273 : _pmetis->options[2] = 15; // random seed (default)
270 327273 : _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 327273 : MeshTools::create_bounding_box(mesh);
283 :
284 327273 : _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 390 : 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 390 : libmesh_assert_equal_to (_pmetis->vtxdist.size(),
297 : cast_int<std::size_t>(mesh.n_processors()+1));
298 390 : libmesh_assert_equal_to (_pmetis->vtxdist[0], 0);
299 :
300 3844200 : for (auto pid : make_range(mesh.n_processors()))
301 : {
302 3518487 : _pmetis->vtxdist[pid+1] = _pmetis->vtxdist[pid] + _n_active_elem_on_proc[pid];
303 3516927 : n_active_elem += _n_active_elem_on_proc[pid];
304 : }
305 390 : 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 780 : std::unordered_map<dof_id_type, dof_id_type> global_index_map;
311 :
312 : {
313 780 : std::vector<dof_id_type> global_index;
314 :
315 : // create the unique mapping for all active elements independent of partitioning
316 : {
317 327663 : MeshBase::const_element_iterator it = mesh.active_elements_begin();
318 327663 : 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 327273 : MeshCommunication().find_global_indices (mesh.comm(),
323 : bbox, it, end,
324 : global_index);
325 :
326 18068385 : for (dof_id_type cnt=0; it != end; ++it)
327 : {
328 17414619 : 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 19794 : libmesh_assert_less (cnt, global_index.size());
333 19794 : libmesh_assert_less (global_index[cnt], n_active_elem);
334 :
335 17434413 : global_index_map.emplace(elem->id(), global_index[cnt++]);
336 : }
337 : }
338 : // really, shouldn't be close!
339 390 : libmesh_assert_less_equal (global_index_map.size(), n_active_elem);
340 390 : 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 327273 : 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 327663 : std::vector<dof_id_type> subdomain_bounds(mesh.n_processors());
355 :
356 328053 : const dof_id_type first_local_elem = _pmetis->vtxdist[mesh.processor_id()];
357 :
358 3844200 : for (auto pid : make_range(mesh.n_processors()))
359 : {
360 780 : dof_id_type tgt_subdomain_size = 0;
361 :
362 : // watch out for the case that n_subdomains < n_processors
363 3516927 : if (pid < static_cast<unsigned int>(_pmetis->nparts))
364 : {
365 2182932 : tgt_subdomain_size = n_active_elem/std::min
366 2182932 : (cast_int<Parmetis::idx_t>(mesh.n_processors()), _pmetis->nparts);
367 :
368 2182152 : if (pid < n_active_elem%_pmetis->nparts)
369 752241 : tgt_subdomain_size++;
370 : }
371 3516927 : if (pid == 0)
372 327273 : subdomain_bounds[0] = tgt_subdomain_size;
373 : else
374 3190434 : subdomain_bounds[pid] = subdomain_bounds[pid-1] + tgt_subdomain_size;
375 : }
376 :
377 390 : libmesh_assert_equal_to (subdomain_bounds.back(), n_active_elem);
378 :
379 : Threads::parallel_for_each
380 327273 : (mesh.active_local_element_stored_range(),
381 3447595 : [first_local_elem, n_active_elem, n_active_local_elem, this,
382 166631 : &global_index_map, &subdomain_bounds](const Elem * elem)
383 : {
384 : const dof_id_type global_index_by_pid =
385 2818277 : libmesh_map_find(_global_index_by_pid_map, elem->id());
386 11834 : libmesh_assert_less (global_index_by_pid, n_active_elem);
387 :
388 2818277 : const dof_id_type local_index =
389 23668 : global_index_by_pid - first_local_elem;
390 :
391 11834 : libmesh_assert_less (local_index, n_active_local_elem);
392 11834 : libmesh_assert_less (local_index, _pmetis->vwgt.size());
393 :
394 11834 : libmesh_ignore(n_active_elem); // unused outside dbg/devel
395 11834 : libmesh_ignore(n_active_local_elem); // unused outside dbg/devel
396 :
397 : // Spline nodes are a special case (storing all the
398 : // unconstrained DoFs in an IGA simulation), but in general
399 : // we'll try to distribute work by expecting it to be roughly
400 : // proportional to DoFs, which are roughly proportional to
401 : // nodes.
402 2818277 : if (elem->type() == NODEELEM &&
403 0 : elem->mapping_type() == RATIONAL_BERNSTEIN_MAP)
404 12789 : _pmetis->vwgt[local_index] = 50;
405 : else
406 2805488 : _pmetis->vwgt[local_index] = elem->n_nodes();
407 :
408 : // find the subdomain this element belongs in
409 : const dof_id_type global_index =
410 2818277 : libmesh_map_find(global_index_map, elem->id());
411 :
412 11834 : libmesh_assert_less (global_index, subdomain_bounds.back());
413 :
414 : const unsigned int subdomain_id =
415 : cast_int<unsigned int>
416 2806443 : (std::distance(subdomain_bounds.begin(),
417 : std::lower_bound(subdomain_bounds.begin(),
418 : subdomain_bounds.end(),
419 11834 : global_index)));
420 11834 : libmesh_assert_less (subdomain_id, _pmetis->nparts);
421 11834 : libmesh_assert_less (local_index, _pmetis->part.size());
422 :
423 2830111 : _pmetis->part[local_index] = subdomain_id;
424 2818277 : });
425 : }
426 327273 : }
427 :
428 :
429 :
430 327273 : void ParmetisPartitioner::build_graph (const MeshBase & mesh)
431 : {
432 1170 : LOG_SCOPE("build_graph()", "ParmetisPartitioner");
433 :
434 : // build the graph in distributed CSR format. Note that
435 : // the edges in the graph will correspond to
436 : // face neighbors
437 390 : const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
438 :
439 327273 : Partitioner::build_graph(mesh);
440 :
441 390 : dof_id_type graph_size_upper_bound = 0;
442 :
443 3145550 : for (auto & row: _dual_graph)
444 2830111 : graph_size_upper_bound += cast_int<dof_id_type>(row.size());
445 :
446 : // Reserve space in the adjacency array
447 327273 : _pmetis->xadj.clear();
448 327273 : _pmetis->xadj.reserve (n_active_local_elem + 1);
449 327273 : _pmetis->adjncy.clear();
450 327273 : _pmetis->adjncy.reserve (graph_size_upper_bound);
451 :
452 : const dof_id_type first_local_elem =
453 328053 : _pmetis->vtxdist[mesh.processor_id()];
454 780 : std::unordered_set<dof_id_type> unique_neighbors;
455 :
456 : // ParMETIS requires a simple graph. Preserve first-occurrence order
457 : // while filtering so valid graph rows remain otherwise unchanged.
458 3145550 : for (auto i : index_range(_dual_graph))
459 : {
460 2830111 : _pmetis->xadj.push_back(cast_int<int>(_pmetis->adjncy.size()));
461 :
462 : const dof_id_type global_index =
463 2818277 : first_local_elem + cast_int<dof_id_type>(i);
464 23668 : const auto & graph_row = _dual_graph[i];
465 11834 : unique_neighbors.clear();
466 :
467 14247737 : for (const auto neighbor : graph_row)
468 11480484 : if (neighbor != global_index &&
469 11429082 : unique_neighbors.insert(neighbor).second)
470 11428018 : _pmetis->adjncy.push_back(neighbor);
471 : }
472 :
473 : // The end of the adjacency array for the last elem
474 327663 : _pmetis->xadj.push_back(cast_int<int>(_pmetis->adjncy.size()));
475 :
476 390 : libmesh_assert_equal_to (_pmetis->xadj.size(), n_active_local_elem+1);
477 390 : libmesh_assert_less_equal (_pmetis->adjncy.size(),
478 : graph_size_upper_bound);
479 327273 : }
480 :
481 : #endif // #ifdef LIBMESH_HAVE_PARMETIS
482 :
483 : } // namespace libMesh
|
