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 : #include "libmesh/partitioner.h"
21 :
22 : // libMesh includes
23 : #include "libmesh/compare_elems_by_level.h"
24 : #include "libmesh/elem.h"
25 : #include "libmesh/enum_to_string.h"
26 : #include "libmesh/int_range.h"
27 : #include "libmesh/libmesh_logging.h"
28 : #include "libmesh/mesh_base.h"
29 : #include "libmesh/mesh_tools.h"
30 : #include "libmesh/mesh_communication.h"
31 : #include "libmesh/parallel_ghost_sync.h"
32 : #include "libmesh/wrapped_petsc.h"
33 : #include "libmesh/boundary_info.h"
34 :
35 : // Subclasses to build()
36 : #include "libmesh/enum_partitioner_type.h"
37 : #include "libmesh/centroid_partitioner.h"
38 : #include "libmesh/hilbert_sfc_partitioner.h"
39 : #include "libmesh/linear_partitioner.h"
40 : #include "libmesh/mapped_subdomain_partitioner.h"
41 : #include "libmesh/metis_partitioner.h"
42 : #include "libmesh/morton_sfc_partitioner.h"
43 : #include "libmesh/parmetis_partitioner.h"
44 : #include "libmesh/sfc_partitioner.h"
45 : #include "libmesh/subdomain_partitioner.h"
46 :
47 : // TIMPI includes
48 : #include "timpi/parallel_implementation.h"
49 : #include "timpi/parallel_sync.h"
50 :
51 : // C/C++ includes
52 : #ifdef LIBMESH_HAVE_PETSC
53 : #include "libmesh/ignore_warnings.h"
54 : #include "petscmat.h"
55 : #include "libmesh/restore_warnings.h"
56 : #include "libmesh/petsc_solver_exception.h"
57 : #endif
58 :
59 :
60 : namespace {
61 :
62 : using namespace libMesh;
63 :
64 : struct CorrectProcIds
65 : {
66 : typedef processor_id_type datum;
67 :
68 144 : CorrectProcIds(MeshBase & _mesh,
69 87375 : std::unordered_set<dof_id_type> & _bad_pids) :
70 87375 : mesh(_mesh), bad_pids(_bad_pids) {}
71 :
72 : MeshBase & mesh;
73 : std::unordered_set<dof_id_type> & bad_pids;
74 :
75 : // ------------------------------------------------------------
76 1547080 : void gather_data (const std::vector<dof_id_type> & ids,
77 : std::vector<datum> & data)
78 : {
79 : // Find the processor id of each requested node
80 1547080 : data.resize(ids.size());
81 :
82 136124319 : for (auto i : index_range(ids))
83 : {
84 : // Look for this point in the mesh and make sure we have a
85 : // good pid for it before sending that on
86 134662944 : if (ids[i] != DofObject::invalid_id &&
87 134557542 : !bad_pids.count(ids[i]))
88 : {
89 134536738 : Node & node = mesh.node_ref(ids[i]);
90 :
91 : // Return the node's correct processor id,
92 134536738 : data[i] = node.processor_id();
93 : }
94 : else
95 40501 : data[i] = DofObject::invalid_processor_id;
96 : }
97 1547080 : }
98 :
99 : // ------------------------------------------------------------
100 1547080 : bool act_on_data (const std::vector<dof_id_type> & ids,
101 : const std::vector<datum> proc_ids)
102 : {
103 462 : bool data_changed = false;
104 :
105 : // Set the ghost node processor ids we've now been informed of
106 136124319 : for (auto i : index_range(ids))
107 : {
108 134577239 : Node & node = mesh.node_ref(ids[i]);
109 :
110 85705 : processor_id_type & proc_id = node.processor_id();
111 134577239 : const processor_id_type new_proc_id = proc_ids[i];
112 :
113 134577239 : if (const auto it = bad_pids.find(ids[i]);
114 134577239 : (new_proc_id != proc_id || it != bad_pids.end()) &&
115 : new_proc_id != DofObject::invalid_processor_id)
116 : {
117 626775 : if (it != bad_pids.end())
118 : {
119 272149 : proc_id = new_proc_id;
120 106 : bad_pids.erase(it);
121 : }
122 : else
123 354626 : proc_id = node.choose_processor_id(proc_id, new_proc_id);
124 :
125 626775 : if (proc_id == new_proc_id)
126 106 : data_changed = true;
127 : }
128 : }
129 :
130 1547080 : return data_changed;
131 : }
132 : };
133 :
134 : }
135 :
136 : namespace libMesh
137 : {
138 :
139 :
140 :
141 : // ------------------------------------------------------------
142 : // Partitioner static data
143 : const dof_id_type Partitioner::communication_blocksize =
144 : dof_id_type(1000000);
145 :
146 :
147 :
148 : // ------------------------------------------------------------
149 : // Partitioner implementation
150 :
151 :
152 0 : PartitionerType Partitioner::type() const
153 : {
154 0 : return INVALID_PARTITIONER;
155 : }
156 :
157 :
158 : std::unique_ptr<Partitioner>
159 418790 : Partitioner::build (const PartitionerType partitioner_type)
160 : {
161 418790 : switch (partitioner_type)
162 : {
163 71 : case CENTROID_PARTITIONER:
164 71 : return std::make_unique<CentroidPartitioner>();
165 213 : case LINEAR_PARTITIONER:
166 213 : return std::make_unique<LinearPartitioner>();
167 0 : case MAPPED_SUBDOMAIN_PARTITIONER:
168 0 : return std::make_unique<MappedSubdomainPartitioner>();
169 69993 : case METIS_PARTITIONER:
170 69993 : return std::make_unique<MetisPartitioner>();
171 348300 : case PARMETIS_PARTITIONER:
172 348300 : return std::make_unique<ParmetisPartitioner>();
173 71 : case HILBERT_SFC_PARTITIONER:
174 71 : return std::make_unique<HilbertSFCPartitioner>();
175 71 : case MORTON_SFC_PARTITIONER:
176 71 : return std::make_unique<MortonSFCPartitioner>();
177 71 : case SFC_PARTITIONER:
178 71 : return std::make_unique<SFCPartitioner>();
179 0 : case SUBDOMAIN_PARTITIONER:
180 0 : return std::make_unique<SubdomainPartitioner>();
181 0 : default:
182 0 : libmesh_error_msg("Invalid partitioner type: " <<
183 : Utility::enum_to_string(partitioner_type));
184 : }
185 : }
186 :
187 :
188 :
189 0 : void Partitioner::partition (MeshBase & mesh)
190 : {
191 0 : this->partition(mesh,mesh.n_processors());
192 0 : }
193 :
194 :
195 :
196 783730 : void Partitioner::partition (MeshBase & mesh,
197 : const unsigned int n)
198 : {
199 15774 : libmesh_parallel_only(mesh.comm());
200 :
201 : // BSK - temporary fix while redistribution is integrated 6/26/2008
202 : // Uncomment this to not repartition in parallel
203 : // if (!mesh.is_serial())
204 : // return;
205 :
206 : // we cannot partition into more pieces than we have
207 : // active elements!
208 : const unsigned int n_parts =
209 : static_cast<unsigned int>
210 783730 : (std::min(mesh.n_active_elem(), static_cast<dof_id_type>(n)));
211 :
212 : // Set the number of partitions in the mesh
213 783730 : mesh.set_n_partitions()=n_parts;
214 :
215 783730 : if (n_parts == 1)
216 : {
217 149914 : this->single_partition (mesh);
218 149914 : return;
219 : }
220 :
221 : // First assign a temporary partitioning to any unpartitioned elements. A spatially
222 : // coherent assignment is only worth its parallel sort for a partitioner that refines what
223 : // it is handed; the other partitioners just need one that is consistent across processors,
224 : // which the cheap traversal-order assignment already is on a replicated mesh.
225 : Partitioner::partition_unpartitioned_elements
226 633816 : (mesh, n_parts,
227 633816 : /*spatial_indexing=*/!(mesh.is_replicated() && this->partitions_from_scratch()));
228 :
229 : // Call the partitioning function
230 633816 : this->_do_partition(mesh,n_parts);
231 :
232 : // Set the parent's processor ids
233 633816 : Partitioner::set_parent_processor_ids(mesh);
234 :
235 : // Redistribute elements if necessary, before setting node processor
236 : // ids, to make sure those will be set consistently
237 633816 : mesh.redistribute();
238 :
239 : #ifdef DEBUG
240 11289 : MeshTools::libmesh_assert_valid_remote_elems(mesh);
241 :
242 : // Messed up elem processor_id()s can leave us without the child
243 : // elements we need to restrict vectors on a distributed mesh
244 11289 : MeshTools::libmesh_assert_valid_procids<Elem>(mesh);
245 : #endif
246 :
247 : // Set the node's processor ids
248 633816 : Partitioner::set_node_processor_ids(mesh);
249 :
250 : #ifdef DEBUG
251 11289 : MeshTools::libmesh_assert_valid_procids<Elem>(mesh);
252 : #endif
253 :
254 : // Give derived Mesh classes a chance to update any cached data to
255 : // reflect the new partitioning
256 633816 : mesh.update_post_partitioning();
257 : }
258 :
259 :
260 :
261 0 : void Partitioner::repartition (MeshBase & mesh)
262 : {
263 0 : this->repartition(mesh,mesh.n_processors());
264 0 : }
265 :
266 :
267 :
268 0 : void Partitioner::repartition (MeshBase & mesh,
269 : const unsigned int n)
270 : {
271 : // we cannot partition into more pieces than we have
272 : // active elements!
273 : const unsigned int n_parts =
274 : static_cast<unsigned int>
275 0 : (std::min(mesh.n_active_elem(), static_cast<dof_id_type>(n)));
276 :
277 : // Set the number of partitions in the mesh
278 0 : mesh.set_n_partitions()=n_parts;
279 :
280 0 : if (n_parts == 1)
281 : {
282 0 : this->single_partition (mesh);
283 0 : return;
284 : }
285 :
286 : // First assign a temporary partitioning to any unpartitioned elements; see partition()
287 : // for why the indexing is chosen this way.
288 : Partitioner::partition_unpartitioned_elements
289 0 : (mesh, n_parts,
290 0 : /*spatial_indexing=*/!(mesh.is_replicated() && this->partitions_from_scratch()));
291 :
292 : // Call the partitioning function
293 0 : this->_do_repartition(mesh,n_parts);
294 :
295 : // Set the parent's processor ids
296 0 : Partitioner::set_parent_processor_ids(mesh);
297 :
298 : // Set the node's processor ids
299 0 : Partitioner::set_node_processor_ids(mesh);
300 : }
301 :
302 :
303 :
304 :
305 :
306 149914 : bool Partitioner::single_partition (MeshBase & mesh)
307 : {
308 : bool changed_pid =
309 295343 : this->single_partition_range(mesh.elements_begin(),
310 440772 : mesh.elements_end());
311 :
312 : // If we have a distributed mesh with an empty rank (or where rank
313 : // 0 has only its own component of a disconnected mesh, I guess),
314 : // that rank might need to be informed of a change.
315 149914 : mesh.comm().max(changed_pid);
316 :
317 : // We may need to redistribute, in case someone (like our unit
318 : // tests) is doing something silly (like moving a whole
319 : // already-distributed mesh back onto rank 0).
320 149914 : if (changed_pid)
321 142062 : mesh.redistribute();
322 :
323 149914 : return changed_pid;
324 : }
325 :
326 :
327 :
328 149914 : bool Partitioner::single_partition_range (MeshBase::element_iterator it,
329 : MeshBase::element_iterator end)
330 : {
331 8970 : LOG_SCOPE("single_partition_range()", "Partitioner");
332 :
333 4485 : bool changed_pid = false;
334 :
335 3219945 : for (auto & elem : as_range(it, end))
336 : {
337 1472380 : if (elem->processor_id())
338 8760 : changed_pid = true;
339 1472380 : elem->processor_id() = 0;
340 :
341 : // Assign all this element's nodes to processor 0 as well.
342 12561548 : for (Node & node : elem->node_ref_range())
343 : {
344 11089168 : if (node.processor_id())
345 40591 : changed_pid = true;
346 11089168 : node.processor_id() = 0;
347 : }
348 140944 : }
349 :
350 154399 : return changed_pid;
351 : }
352 :
353 0 : void Partitioner::partition_unpartitioned_elements (MeshBase & mesh)
354 : {
355 0 : Partitioner::partition_unpartitioned_elements(mesh, mesh.n_processors());
356 0 : }
357 :
358 :
359 :
360 0 : void Partitioner::partition_unpartitioned_elements (MeshBase & mesh,
361 : const unsigned int n_subdomains)
362 : {
363 0 : Partitioner::partition_unpartitioned_elements(mesh, n_subdomains,
364 : /*spatial_indexing=*/true);
365 0 : }
366 :
367 :
368 :
369 633816 : void Partitioner::partition_unpartitioned_elements (MeshBase & mesh,
370 : const unsigned int n_subdomains,
371 : const bool spatial_indexing)
372 : {
373 633816 : MeshBase::element_iterator it = mesh.unpartitioned_elements_begin();
374 633816 : const MeshBase::element_iterator end = mesh.unpartitioned_elements_end();
375 :
376 1256343 : const dof_id_type n_unpartitioned_elements = MeshTools::n_elem (it, end);
377 :
378 : // the unpartitioned elements must exist on all processors. If the range is empty on one
379 : // it is empty on all, and we can quit right here.
380 633816 : if (!n_unpartitioned_elements)
381 3679 : return;
382 :
383 : // find the target subdomain sizes
384 276043 : std::vector<dof_id_type> subdomain_bounds(mesh.n_processors());
385 :
386 2871756 : for (auto pid : make_range(mesh.n_processors()))
387 : {
388 15220 : dof_id_type tgt_subdomain_size = 0;
389 :
390 : // watch out for the case that n_subdomains < n_processors
391 2610933 : if (pid < n_subdomains)
392 : {
393 1454770 : tgt_subdomain_size = n_unpartitioned_elements/n_subdomains;
394 :
395 1454770 : if (pid < n_unpartitioned_elements%n_subdomains)
396 420410 : tgt_subdomain_size++;
397 :
398 : }
399 :
400 : //libMesh::out << "pid, #= " << pid << ", " << tgt_subdomain_size << std::endl;
401 2610933 : if (pid == 0)
402 260823 : subdomain_bounds[0] = tgt_subdomain_size;
403 : else
404 2365330 : subdomain_bounds[pid] = subdomain_bounds[pid-1] + tgt_subdomain_size;
405 : }
406 :
407 7610 : libmesh_assert_equal_to (subdomain_bounds.back(), n_unpartitioned_elements);
408 :
409 : // create the unique mapping for all unpartitioned elements independent of partitioning
410 : // determine the global indexing for all the unpartitioned elements
411 15220 : std::vector<dof_id_type> global_indices;
412 :
413 : // Calling this on all processors a unique range in [0,n_unpartitioned_elements) is constructed.
414 : // Only the indices for the elements we pass in are returned in the array.
415 : // Without the spatial indexing an element's position in the traversal of the range is used
416 : // instead, which a replicated mesh reaches in the same order on every processor and which
417 : // therefore indexes [0, n_unpartitioned_elements) just as uniquely, at no cost.
418 260823 : if (spatial_indexing)
419 223545 : MeshCommunication().find_global_indices (mesh.comm(),
420 447090 : MeshTools::create_bounding_box(mesh), it, end,
421 : global_indices);
422 : else
423 7398 : libmesh_assert (mesh.is_replicated());
424 :
425 7610 : dof_id_type cnt=0;
426 14001835 : for (auto & elem : as_range(it, end))
427 : {
428 374320 : libmesh_assert (!spatial_indexing || cnt < global_indices.size());
429 : const dof_id_type global_index =
430 13487799 : spatial_indexing ? global_indices[cnt++] : cnt++;
431 :
432 374320 : libmesh_assert_less (global_index, subdomain_bounds.back());
433 374320 : libmesh_assert_less (global_index, n_unpartitioned_elements);
434 :
435 : const processor_id_type subdomain_id =
436 : cast_int<processor_id_type>
437 13113479 : (std::distance(subdomain_bounds.begin(),
438 : std::upper_bound(subdomain_bounds.begin(),
439 : subdomain_bounds.end(),
440 374320 : global_index)));
441 374320 : libmesh_assert_less (subdomain_id, n_subdomains);
442 :
443 13487799 : elem->processor_id() = subdomain_id;
444 : //libMesh::out << "assigning " << global_index << " to " << subdomain_id << std::endl;
445 245603 : }
446 : }
447 :
448 :
449 :
450 633816 : void Partitioner::set_parent_processor_ids(MeshBase & mesh)
451 : {
452 : // Ignore the parameter when !LIBMESH_ENABLE_AMR
453 11289 : libmesh_ignore(mesh);
454 :
455 33867 : LOG_SCOPE("set_parent_processor_ids()", "Partitioner");
456 :
457 : #ifdef LIBMESH_ENABLE_AMR
458 :
459 : // If the mesh is serial we have access to all the elements,
460 : // in particular all the active ones. We can therefore set
461 : // the parent processor ids indirectly through their children, and
462 : // set the subactive processor ids while examining their active
463 : // ancestors.
464 : // By convention a parent is assigned to the minimum processor
465 : // of all its children, and a subactive is assigned to the processor
466 : // of its active ancestor.
467 633816 : if (mesh.is_serial())
468 : {
469 52786987 : for (auto & elem : mesh.active_element_ptr_range())
470 : {
471 : // First set descendents
472 2268744 : std::vector<Elem *> subactive_family;
473 26975115 : elem->total_family_tree(subactive_family);
474 55210266 : for (const auto & f : subactive_family)
475 28235151 : f->processor_id() = elem->processor_id();
476 :
477 : // Then set ancestors
478 26975115 : Elem * parent = elem->parent();
479 :
480 60375333 : while (parent)
481 : {
482 : // invalidate the parent id, otherwise the min below
483 : // will not work if the current parent id is less
484 : // than all the children!
485 1988368 : parent->invalidate_processor_id();
486 :
487 173394452 : for (auto & child : parent->child_ref_range())
488 : {
489 9496964 : libmesh_assert(!child.is_remote());
490 9496964 : libmesh_assert_not_equal_to (child.processor_id(), DofObject::invalid_processor_id);
491 139994234 : parent->processor_id() = std::min(parent->processor_id(),
492 9496964 : child.processor_id());
493 : }
494 3976736 : parent = parent->parent();
495 : }
496 536012 : }
497 : }
498 :
499 : // When the mesh is parallel we cannot guarantee that parents have access to
500 : // all their children.
501 : else
502 : {
503 : // Setting subactive processor ids is easy: we can guarantee
504 : // that children have access to all their parents.
505 :
506 : // Loop over all the active elements in the mesh
507 10877280 : for (auto & child : mesh.active_element_ptr_range())
508 : {
509 22684 : std::vector<Elem *> subactive_family;
510 5374561 : child->total_family_tree(subactive_family);
511 10752863 : for (const auto & f : subactive_family)
512 5378302 : f->processor_id() = child->processor_id();
513 75226 : }
514 :
515 : // When the mesh is parallel we cannot guarantee that parents have access to
516 : // all their children.
517 :
518 : // We will use a brute-force approach here. Each processor finds its parent
519 : // elements and sets the parent pid to the minimum of its
520 : // semilocal descendants.
521 : // A global reduction is then performed to make sure the true minimum is found.
522 : // As noted, this is required because we cannot guarantee that a parent has
523 : // access to all its children on any single processor.
524 130 : libmesh_parallel_only(mesh.comm());
525 130 : libmesh_assert(MeshTools::n_elem(mesh.unpartitioned_elements_begin(),
526 : mesh.unpartitioned_elements_end()) == 0);
527 :
528 75486 : const dof_id_type max_elem_id = mesh.max_elem_id();
529 :
530 : std::vector<processor_id_type>
531 260 : parent_processor_ids (std::min(communication_blocksize,
532 75616 : max_elem_id));
533 :
534 150972 : for (dof_id_type blk=0, last_elem_id=0; last_elem_id<max_elem_id; blk++)
535 : {
536 75486 : last_elem_id =
537 75616 : std::min(static_cast<dof_id_type>((blk+1)*communication_blocksize),
538 130 : max_elem_id);
539 75486 : const dof_id_type first_elem_id = blk*communication_blocksize;
540 :
541 130 : std::fill (parent_processor_ids.begin(),
542 : parent_processor_ids.end(),
543 : DofObject::invalid_processor_id);
544 :
545 : // first build up local contributions to parent_processor_ids
546 130 : bool have_parent_in_block = false;
547 :
548 151834 : for (auto & parent : as_range(mesh.ancestor_elements_begin(),
549 1674418 : mesh.ancestor_elements_end()))
550 : {
551 686863 : const dof_id_type parent_idx = parent->id();
552 992 : libmesh_assert_less (parent_idx, max_elem_id);
553 :
554 687855 : if ((parent_idx >= first_elem_id) &&
555 686863 : (parent_idx < last_elem_id))
556 : {
557 992 : have_parent_in_block = true;
558 686863 : processor_id_type parent_pid = DofObject::invalid_processor_id;
559 :
560 992 : std::vector<const Elem *> active_family;
561 686863 : parent->active_family_tree(active_family);
562 7909677 : for (const auto & f : active_family)
563 8010075 : parent_pid = std::min (parent_pid, f->processor_id());
564 :
565 686863 : const dof_id_type packed_idx = parent_idx - first_elem_id;
566 992 : libmesh_assert_less (packed_idx, parent_processor_ids.size());
567 :
568 687855 : parent_processor_ids[packed_idx] = parent_pid;
569 : }
570 75226 : }
571 :
572 : // then find the global minimum
573 75486 : mesh.comm().min (parent_processor_ids);
574 :
575 : // and assign the ids, if we have a parent in this block.
576 75486 : if (have_parent_in_block)
577 48956 : for (auto & parent : as_range(mesh.ancestor_elements_begin(),
578 1468662 : mesh.ancestor_elements_end()))
579 : {
580 686863 : const dof_id_type parent_idx = parent->id();
581 :
582 687855 : if ((parent_idx >= first_elem_id) &&
583 685871 : (parent_idx < last_elem_id))
584 : {
585 686863 : const dof_id_type packed_idx = parent_idx - first_elem_id;
586 992 : libmesh_assert_less (packed_idx, parent_processor_ids.size());
587 :
588 : const processor_id_type parent_pid =
589 686863 : parent_processor_ids[packed_idx];
590 :
591 992 : libmesh_assert_not_equal_to (parent_pid, DofObject::invalid_processor_id);
592 :
593 686863 : parent->processor_id() = parent_pid;
594 : }
595 23904 : }
596 : }
597 : }
598 :
599 : #endif // LIBMESH_ENABLE_AMR
600 633816 : }
601 :
602 : void
603 0 : Partitioner::processor_pairs_to_interface_nodes(MeshBase & mesh,
604 : std::map<std::pair<processor_id_type, processor_id_type>, std::set<dof_id_type>> & processor_pair_to_nodes)
605 : {
606 : // This function must be run on all processors at once
607 0 : libmesh_parallel_only(mesh.comm());
608 :
609 0 : processor_pair_to_nodes.clear();
610 :
611 0 : std::set<dof_id_type> mynodes;
612 0 : std::set<dof_id_type> neighbor_nodes;
613 0 : std::vector<dof_id_type> common_nodes;
614 :
615 : // Loop over all the active elements
616 0 : for (auto & elem : mesh.active_element_ptr_range())
617 : {
618 0 : libmesh_assert(elem);
619 :
620 0 : libmesh_assert_not_equal_to (elem->processor_id(), DofObject::invalid_processor_id);
621 :
622 0 : auto n_nodes = elem->n_nodes();
623 :
624 : // prepare data for this element
625 0 : mynodes.clear();
626 0 : neighbor_nodes.clear();
627 0 : common_nodes.clear();
628 :
629 0 : for (unsigned int inode = 0; inode < n_nodes; inode++)
630 0 : mynodes.insert(elem->node_id(inode));
631 :
632 0 : for (auto i : elem->side_index_range())
633 : {
634 0 : auto neigh = elem->neighbor_ptr(i);
635 0 : if (neigh && !neigh->is_remote() && neigh->processor_id() != elem->processor_id())
636 : {
637 0 : neighbor_nodes.clear();
638 0 : common_nodes.clear();
639 0 : auto neigh_n_nodes = neigh->n_nodes();
640 0 : for (unsigned int inode = 0; inode < neigh_n_nodes; inode++)
641 0 : neighbor_nodes.insert(neigh->node_id(inode));
642 :
643 0 : std::set_intersection(mynodes.begin(), mynodes.end(),
644 : neighbor_nodes.begin(), neighbor_nodes.end(),
645 0 : std::back_inserter(common_nodes));
646 :
647 0 : auto & map_set = processor_pair_to_nodes[std::make_pair(std::min(elem->processor_id(), neigh->processor_id()),
648 0 : std::max(elem->processor_id(), neigh->processor_id()))];
649 0 : for (auto global_node_id : common_nodes)
650 0 : map_set.insert(global_node_id);
651 : }
652 : }
653 0 : }
654 0 : }
655 :
656 0 : void Partitioner::set_interface_node_processor_ids_linear(MeshBase & mesh)
657 : {
658 : // This function must be run on all processors at once
659 0 : libmesh_parallel_only(mesh.comm());
660 :
661 0 : std::map<std::pair<processor_id_type, processor_id_type>, std::set<dof_id_type>> processor_pair_to_nodes;
662 :
663 0 : processor_pairs_to_interface_nodes(mesh, processor_pair_to_nodes);
664 :
665 0 : for (auto & pmap : processor_pair_to_nodes)
666 : {
667 0 : std::size_t n_own_nodes = pmap.second.size()/2, i = 0;
668 :
669 0 : for (dof_id_type id : pmap.second)
670 : {
671 0 : auto & node = mesh.node_ref(id);
672 0 : if (i <= n_own_nodes)
673 0 : node.processor_id() = pmap.first.first;
674 : else
675 0 : node.processor_id() = pmap.first.second;
676 0 : i++;
677 : }
678 : }
679 0 : }
680 :
681 0 : void Partitioner::set_interface_node_processor_ids_BFS(MeshBase & mesh)
682 : {
683 : // This function must be run on all processors at once
684 0 : libmesh_parallel_only(mesh.comm());
685 :
686 : // I see occasional consistency failures when using this on a
687 : // distributed mesh
688 : libmesh_experimental();
689 :
690 0 : std::map<std::pair<processor_id_type, processor_id_type>, std::set<dof_id_type>> processor_pair_to_nodes;
691 :
692 0 : processor_pairs_to_interface_nodes(mesh, processor_pair_to_nodes);
693 :
694 0 : std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
695 :
696 0 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
697 :
698 0 : std::vector<const Node *> neighbors;
699 0 : std::set<dof_id_type> neighbors_order;
700 0 : std::vector<dof_id_type> common_nodes;
701 0 : std::queue<dof_id_type> nodes_queue;
702 0 : std::set<dof_id_type> visted_nodes;
703 :
704 0 : for (auto & pmap : processor_pair_to_nodes)
705 : {
706 0 : std::size_t n_own_nodes = pmap.second.size()/2;
707 :
708 : // Initialize node assignment
709 0 : for (dof_id_type id : pmap.second)
710 0 : mesh.node_ref(id).processor_id() = pmap.first.second;
711 :
712 0 : visted_nodes.clear();
713 0 : for (dof_id_type id : pmap.second)
714 : {
715 0 : mesh.node_ref(id).processor_id() = pmap.first.second;
716 :
717 0 : if (visted_nodes.count(id))
718 0 : continue;
719 : else
720 : {
721 0 : nodes_queue.push(id);
722 0 : visted_nodes.insert(id);
723 0 : if (visted_nodes.size() >= n_own_nodes)
724 0 : break;
725 : }
726 :
727 0 : while (!nodes_queue.empty())
728 : {
729 0 : auto & node = mesh.node_ref(nodes_queue.front());
730 0 : nodes_queue.pop();
731 :
732 0 : neighbors.clear();
733 0 : MeshTools::find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
734 0 : neighbors_order.clear();
735 0 : for (auto & neighbor : neighbors)
736 0 : neighbors_order.insert(neighbor->id());
737 :
738 0 : common_nodes.clear();
739 0 : std::set_intersection(pmap.second.begin(), pmap.second.end(),
740 : neighbors_order.begin(), neighbors_order.end(),
741 0 : std::back_inserter(common_nodes));
742 :
743 0 : for (auto c_node : common_nodes)
744 0 : if (!visted_nodes.count(c_node))
745 : {
746 0 : nodes_queue.push(c_node);
747 0 : visted_nodes.insert(c_node);
748 0 : if (visted_nodes.size() >= n_own_nodes)
749 0 : goto queue_done;
750 : }
751 :
752 0 : if (visted_nodes.size() >= n_own_nodes)
753 0 : goto queue_done;
754 : }
755 : }
756 0 : queue_done:
757 0 : for (auto node : visted_nodes)
758 0 : mesh.node_ref(node).processor_id() = pmap.first.first;
759 : }
760 0 : }
761 :
762 0 : void Partitioner::set_interface_node_processor_ids_petscpartitioner(MeshBase & mesh)
763 : {
764 0 : libmesh_ignore(mesh); // Only used if LIBMESH_HAVE_PETSC
765 :
766 : // This function must be run on all processors at once
767 0 : libmesh_parallel_only(mesh.comm());
768 :
769 : #ifdef LIBMESH_HAVE_PETSC
770 0 : std::map<std::pair<processor_id_type, processor_id_type>, std::set<dof_id_type>> processor_pair_to_nodes;
771 :
772 0 : processor_pairs_to_interface_nodes(mesh, processor_pair_to_nodes);
773 :
774 0 : std::vector<std::vector<const Elem *>> nodes_to_elem_map;
775 :
776 0 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
777 :
778 0 : std::vector<const Node *> neighbors;
779 0 : std::set<dof_id_type> neighbors_order;
780 0 : std::vector<dof_id_type> common_nodes;
781 :
782 0 : std::vector<dof_id_type> rows;
783 0 : std::vector<dof_id_type> cols;
784 :
785 0 : std::map<dof_id_type, dof_id_type> global_to_local;
786 :
787 0 : for (auto & pmap : processor_pair_to_nodes)
788 : {
789 0 : unsigned int i = 0;
790 :
791 0 : rows.clear();
792 0 : rows.resize(pmap.second.size()+1);
793 0 : cols.clear();
794 0 : for (dof_id_type id : pmap.second)
795 0 : global_to_local[id] = i++;
796 :
797 0 : i = 0;
798 0 : for (auto id : pmap.second)
799 : {
800 0 : auto & node = mesh.node_ref(id);
801 0 : neighbors.clear();
802 0 : MeshTools::find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
803 0 : neighbors_order.clear();
804 0 : for (auto & neighbor : neighbors)
805 0 : neighbors_order.insert(neighbor->id());
806 :
807 0 : common_nodes.clear();
808 0 : std::set_intersection(pmap.second.begin(), pmap.second.end(),
809 : neighbors_order.begin(), neighbors_order.end(),
810 0 : std::back_inserter(common_nodes));
811 :
812 0 : rows[i+1] = rows[i] + cast_int<dof_id_type>(common_nodes.size());
813 :
814 0 : for (auto c_node : common_nodes)
815 0 : cols.push_back(global_to_local[c_node]);
816 :
817 0 : i++;
818 : }
819 :
820 : // Next we construct an IS from a MatPartitioning
821 0 : WrappedPetsc<IS> is;
822 : {
823 : PetscInt *adj_i, *adj_j;
824 0 : LibmeshPetscCall2(mesh.comm(), PetscCalloc1(rows.size(), &adj_i));
825 0 : LibmeshPetscCall2(mesh.comm(), PetscCalloc1(cols.size(), &adj_j));
826 0 : PetscInt rows_size = cast_int<PetscInt>(rows.size());
827 0 : for (PetscInt ii=0; ii<rows_size; ii++)
828 0 : adj_i[ii] = rows[ii];
829 :
830 0 : PetscInt cols_size = cast_int<PetscInt>(cols.size());
831 0 : for (PetscInt ii=0; ii<cols_size; ii++)
832 0 : adj_j[ii] = cols[ii];
833 :
834 0 : const PetscInt sz = cast_int<PetscInt>(pmap.second.size());
835 :
836 : // Create sparse matrix representing an adjacency list
837 0 : WrappedPetsc<Mat> adj;
838 0 : LibmeshPetscCall2(mesh.comm(), MatCreateMPIAdj(PETSC_COMM_SELF, sz, sz, adj_i, adj_j, nullptr, adj.get()));
839 :
840 : // Create MatPartitioning object
841 0 : WrappedPetsc<MatPartitioning> part;
842 0 : LibmeshPetscCall2(mesh.comm(), MatPartitioningCreate(PETSC_COMM_SELF, part.get()));
843 :
844 : // Apply MatPartitioning, storing results in "is"
845 0 : LibmeshPetscCall2(mesh.comm(), MatPartitioningSetAdjacency(part, adj));
846 0 : LibmeshPetscCall2(mesh.comm(), MatPartitioningSetNParts(part, 2));
847 0 : LibmeshPetscCall2(mesh.comm(), PetscObjectSetOptionsPrefix((PetscObject)(*part), "balance_"));
848 0 : LibmeshPetscCall2(mesh.comm(), MatPartitioningSetFromOptions(part));
849 0 : LibmeshPetscCall2(mesh.comm(), MatPartitioningApply(part, is.get()));
850 : }
851 :
852 : PetscInt local_size;
853 : const PetscInt *indices;
854 0 : LibmeshPetscCall2(mesh.comm(), ISGetLocalSize(is, &local_size));
855 0 : LibmeshPetscCall2(mesh.comm(), ISGetIndices(is, &indices));
856 :
857 0 : i = 0;
858 0 : for (auto id : pmap.second)
859 : {
860 0 : auto & node = mesh.node_ref(id);
861 0 : if (indices[i])
862 0 : node.processor_id() = pmap.first.second;
863 : else
864 0 : node.processor_id() = pmap.first.first;
865 :
866 0 : i++;
867 : }
868 0 : LibmeshPetscCall2(mesh.comm(), ISRestoreIndices(is, &indices));
869 : }
870 : #else
871 0 : libmesh_error_msg("PETSc is required");
872 : #endif
873 0 : }
874 :
875 :
876 672385 : void Partitioner::set_node_processor_ids(MeshBase & mesh)
877 : {
878 39771 : LOG_SCOPE("set_node_processor_ids()","Partitioner");
879 :
880 : // This function must be run on all processors at once
881 13257 : libmesh_parallel_only(mesh.comm());
882 :
883 : // If we have any unpartitioned elements at this
884 : // stage there is a problem
885 13257 : libmesh_assert (MeshTools::n_elem(mesh.unpartitioned_elements_begin(),
886 : mesh.unpartitioned_elements_end()) == 0);
887 :
888 : // Start from scratch here: nodes we used to own may not be
889 : // eligible for us to own any more.
890 98724729 : for (auto & node : mesh.node_ptr_range())
891 : {
892 97393216 : node->processor_id() = DofObject::invalid_processor_id;
893 645871 : }
894 :
895 : // Loop over all the active elements
896 79346111 : for (auto & elem : mesh.active_element_ptr_range())
897 : {
898 1779932 : libmesh_assert(elem);
899 :
900 1779932 : libmesh_assert_not_equal_to (elem->processor_id(), DofObject::invalid_processor_id);
901 :
902 : // Consider updating the processor id on this element's nodes
903 356495846 : for (Node & node : elem->node_ref_range())
904 : {
905 12667995 : processor_id_type & pid = node.processor_id();
906 313928683 : pid = node.choose_processor_id(pid, elem->processor_id());
907 : }
908 645871 : }
909 :
910 : // How we finish off the node partitioning depends on our command
911 : // line options.
912 :
913 : const bool load_balanced_nodes_linear =
914 672385 : libMesh::on_command_line ("--load-balanced-nodes-linear");
915 :
916 : const bool load_balanced_nodes_bfs =
917 672385 : libMesh::on_command_line ("--load-balanced-nodes-bfs");
918 :
919 : const bool load_balanced_nodes_petscpartition =
920 672385 : libMesh::on_command_line ("--load-balanced-nodes-petscpartitioner");
921 :
922 672385 : unsigned int n_load_balance_options = load_balanced_nodes_linear;
923 672385 : n_load_balance_options += load_balanced_nodes_bfs;
924 672385 : n_load_balance_options += load_balanced_nodes_petscpartition;
925 672385 : libmesh_error_msg_if(n_load_balance_options > 1,
926 : "Cannot perform more than one load balancing type at a time");
927 :
928 672385 : if (load_balanced_nodes_linear)
929 0 : set_interface_node_processor_ids_linear(mesh);
930 672385 : else if (load_balanced_nodes_bfs)
931 0 : set_interface_node_processor_ids_BFS(mesh);
932 672385 : else if (load_balanced_nodes_petscpartition)
933 0 : set_interface_node_processor_ids_petscpartitioner(mesh);
934 :
935 : // Node balancing algorithm will response to assign owned nodes.
936 : // We still need to sync PIDs
937 : {
938 : // For inactive elements, we will have already gotten most of
939 : // these nodes, *except* for the case of a parent with a subset
940 : // of active descendants which are remote elements. In that
941 : // case some of the parent nodes will not have been properly
942 : // handled yet on our processor.
943 : //
944 : // We don't want to inadvertently give one of them an incorrect
945 : // processor id, but if we're not in serial then we have to
946 : // assign them temporary pids to make querying work, so we'll
947 : // save our *valid* pids before assigning temporaries.
948 : //
949 : // Even in serial we'll want to check and make sure we're not
950 : // overwriting valid active node pids with pids from subactive
951 : // elements.
952 26514 : std::unordered_set<dof_id_type> bad_pids;
953 :
954 189037865 : for (auto & node : mesh.node_ptr_range())
955 97393216 : if (node->processor_id() == DofObject::invalid_processor_id)
956 2648545 : bad_pids.insert(node->id());
957 :
958 : // If we assign our temporary ids by looping from finer elements
959 : // to coarser elements, we'll always get an id from the finest
960 : // ghost element we can see, which will usually be "closer" to
961 : // the true processor we want to query and so will reduce query
962 : // cycles that don't reach that processor.
963 :
964 : // But we can still end up with a query cycle that dead-ends, so
965 : // we need to prepare a "push" communication step here.
966 :
967 672385 : const bool is_serial = mesh.is_serial();
968 : std::unordered_map
969 : <processor_id_type,
970 : std::unordered_map<dof_id_type, processor_id_type>>
971 26514 : potential_pids;
972 :
973 672385 : const unsigned int n_levels = MeshTools::n_levels(mesh);
974 1566314 : for (unsigned int level = n_levels; level > 0; --level)
975 : {
976 3935399 : for (auto & elem : as_range(mesh.level_elements_begin(level-1),
977 98201016 : mesh.level_elements_end(level-1)))
978 : {
979 2169376 : libmesh_assert_not_equal_to (elem->processor_id(),
980 : DofObject::invalid_processor_id);
981 :
982 48419173 : const processor_id_type elem_pid = elem->processor_id();
983 :
984 : // Consider updating the processor id on this element's nodes
985 398205429 : for (Node & node : elem->node_ref_range())
986 : {
987 14612985 : processor_id_type & pid = node.processor_id();
988 349786256 : if (bad_pids.count(node.id()))
989 5398401 : pid = node.choose_processor_id(pid, elem_pid);
990 344387855 : else if (!is_serial)
991 61187363 : potential_pids[elem_pid][node.id()] = pid;
992 : }
993 850259 : }
994 : }
995 :
996 672385 : if (!is_serial)
997 : {
998 : std::unordered_map
999 : <processor_id_type,
1000 : std::vector<std::pair<dof_id_type, processor_id_type>>>
1001 144 : potential_pids_vecs;
1002 :
1003 643797 : for (auto & pair : potential_pids)
1004 556422 : potential_pids_vecs[pair.first].assign(pair.second.begin(), pair.second.end());
1005 :
1006 : auto pids_action_functor =
1007 555886 : [& mesh, & bad_pids]
1008 : (processor_id_type /* src_pid */,
1009 23513389 : const std::vector<std::pair<dof_id_type, processor_id_type>> & data)
1010 : {
1011 23964250 : for (auto pair : data)
1012 : {
1013 23407828 : Node & node = mesh.node_ref(pair.first);
1014 35089 : processor_id_type & pid = node.processor_id();
1015 23407828 : if (const auto it = bad_pids.find(pair.first);
1016 35089 : it != bad_pids.end())
1017 : {
1018 15678 : pid = pair.second;
1019 26 : bad_pids.erase(it);
1020 : }
1021 : else
1022 23392150 : pid = node.choose_processor_id(pid, pair.second);
1023 : }
1024 87767 : };
1025 :
1026 : Parallel::push_parallel_vector_data
1027 87375 : (mesh.comm(), potential_pids_vecs, pids_action_functor);
1028 :
1029 : // Using default libMesh options, we'll just need to sync
1030 : // between processors now. The catch here is that we can't
1031 : // initially trust Node::choose_processor_id() because some
1032 : // of those node processor ids are the temporary ones.
1033 144 : CorrectProcIds correct_pids(mesh, bad_pids);
1034 : Parallel::sync_node_data_by_element_id
1035 261837 : (mesh, mesh.elements_begin(), mesh.elements_end(),
1036 87231 : Parallel::SyncEverything(), Parallel::SyncEverything(),
1037 : correct_pids);
1038 :
1039 : // But once we've got all the non-temporary pids synced, we
1040 : // may need to sync again to get any pids on nodes only
1041 : // connected to subactive elements, for which *only*
1042 : // "temporary" pids are possible.
1043 144 : bad_pids.clear();
1044 : Parallel::sync_node_data_by_element_id
1045 174750 : (mesh,
1046 262269 : mesh.elements_begin(), mesh.elements_end(),
1047 87375 : Parallel::SyncEverything(), Parallel::SyncEverything(),
1048 : correct_pids);
1049 : }
1050 : }
1051 :
1052 : // We can't assert that all nodes are connected to elements, because
1053 : // a DistributedMesh with NodeConstraints might have pulled in some
1054 : // remote nodes solely for evaluating those constraints.
1055 : // MeshTools::libmesh_assert_connected_nodes(mesh);
1056 :
1057 : #ifdef DEBUG
1058 13257 : MeshTools::libmesh_assert_valid_procids<Node>(mesh);
1059 : //MeshTools::libmesh_assert_canonical_node_procids(mesh);
1060 : #endif
1061 672385 : }
1062 :
1063 :
1064 :
1065 : struct SyncLocalIDs
1066 : {
1067 : typedef dof_id_type datum;
1068 :
1069 : typedef std::unordered_map<dof_id_type, dof_id_type> map_type;
1070 :
1071 624762 : SyncLocalIDs(map_type & _id_map) : id_map(_id_map) {}
1072 :
1073 : map_type & id_map;
1074 :
1075 3230402 : void gather_data (const std::vector<dof_id_type> & ids,
1076 : std::vector<datum> & local_ids) const
1077 : {
1078 3231152 : local_ids.resize(ids.size());
1079 :
1080 32181282 : for (auto i : index_range(ids))
1081 28950880 : local_ids[i] = id_map[ids[i]];
1082 3230402 : }
1083 :
1084 3230402 : void act_on_data (const std::vector<dof_id_type> & ids,
1085 : const std::vector<datum> & local_ids)
1086 : {
1087 32181282 : for (auto i : index_range(local_ids))
1088 28966784 : id_map[ids[i]] = local_ids[i];
1089 3230402 : }
1090 : };
1091 :
1092 624762 : void Partitioner::_find_global_index_by_pid_map(const MeshBase & mesh)
1093 : {
1094 624762 : const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
1095 :
1096 : // Find the number of active elements on each processor. We cannot use
1097 : // mesh.n_active_elem_on_proc(pid) since that only returns the number of
1098 : // elements assigned to pid which are currently stored on the calling
1099 : // processor. This will not in general be correct for parallel meshes
1100 : // when (pid!=mesh.processor_id()).
1101 1528 : auto n_proc = mesh.n_processors();
1102 624762 : _n_active_elem_on_proc.resize(n_proc);
1103 624762 : mesh.comm().allgather(n_active_local_elem, _n_active_elem_on_proc);
1104 :
1105 626290 : std::vector<dof_id_type> n_active_elem_before_proc(mesh.n_processors());
1106 :
1107 6713270 : for (auto i : make_range(n_proc-1))
1108 6088508 : n_active_elem_before_proc[i+1] =
1109 6090036 : n_active_elem_before_proc[i] + _n_active_elem_on_proc[i];
1110 :
1111 : libMesh::BoundingBox bbox =
1112 624762 : MeshTools::create_bounding_box(mesh);
1113 :
1114 764 : _global_index_by_pid_map.clear();
1115 :
1116 : // create the mapping which is contiguous by processor
1117 624762 : MeshCommunication().find_local_indices (bbox,
1118 1249524 : mesh.active_local_elements_begin(),
1119 1248760 : mesh.active_local_elements_end(),
1120 624762 : _global_index_by_pid_map);
1121 :
1122 764 : SyncLocalIDs sync(_global_index_by_pid_map);
1123 :
1124 : Parallel::sync_dofobject_data_by_id
1125 1248760 : (mesh.comm(), mesh.active_elements_begin(), mesh.active_elements_end(), sync);
1126 :
1127 35804574 : for (const auto & elem : mesh.active_element_ptr_range())
1128 : {
1129 34555814 : const processor_id_type pid = elem->processor_id();
1130 39556 : libmesh_assert_less (_global_index_by_pid_map[elem->id()], _n_active_elem_on_proc[pid]);
1131 :
1132 34595370 : _global_index_by_pid_map[elem->id()] += n_active_elem_before_proc[pid];
1133 623234 : }
1134 624762 : }
1135 :
1136 312381 : void Partitioner::build_graph (const MeshBase & mesh)
1137 : {
1138 1146 : LOG_SCOPE("build_graph()", "Partitioner");
1139 :
1140 382 : const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
1141 :
1142 : // If we have boundary elements in this mesh, we want to account for
1143 : // the connectivity between them and interior elements. We can find
1144 : // interior elements from boundary elements, but we need to build up
1145 : // a lookup map to do the reverse.
1146 : typedef std::unordered_multimap<const Elem *, const Elem *> map_type;
1147 764 : map_type interior_to_boundary_map;
1148 :
1149 : // If we have spline nodes in this mesh, we want to account for the
1150 : // connectivity between them and integration elements. We can find
1151 : // spline nodes from integration elements, but need a reverse map
1152 : // {integration_elements} = elems_constrained_by[spline_nodeelem]
1153 764 : map_type elems_constrained_by;
1154 :
1155 382 : const auto & mesh_constrained_nodes = mesh.get_constraint_rows();
1156 :
1157 52398767 : for (const Elem * elem : mesh.active_element_ptr_range())
1158 : {
1159 17277907 : if (!mesh_constrained_nodes.empty()) // quick test for non-IGA cases
1160 : {
1161 0 : const auto end_it = mesh_constrained_nodes.end();
1162 :
1163 : // Use a set to avoid duplicates. Use a well-defined method
1164 : // of ordering that set to make debugging easier.
1165 0 : std::set<const Elem *, CompareElemIdsByLevel> constraining_elems;
1166 1209187 : for (const Node & node : elem->node_ref_range())
1167 : {
1168 980466 : if (const auto row_it = mesh_constrained_nodes.find(&node);
1169 0 : row_it != end_it)
1170 3946461 : for (const auto & [pr, coef] : row_it->second)
1171 : {
1172 0 : libmesh_ignore(coef); // avoid gcc 7 warning
1173 3092771 : constraining_elems.insert(pr.first);
1174 : }
1175 : }
1176 1583959 : for (const Elem * constraining_elem : constraining_elems)
1177 0 : elems_constrained_by.emplace(constraining_elem, elem);
1178 : }
1179 :
1180 : // If we don't have an interior_parent and we don't have any
1181 : // constrained nodes then there's nothing else to look up.
1182 17277907 : if (elem->interior_parent())
1183 : {
1184 : // get all relevant interior elements
1185 0 : std::set<const Elem *> neighbor_set;
1186 1094 : elem->find_interior_neighbors(neighbor_set);
1187 :
1188 2188 : for (const auto & neighbor : neighbor_set)
1189 0 : interior_to_boundary_map.emplace(neighbor, elem);
1190 : }
1191 311617 : }
1192 :
1193 : #ifdef LIBMESH_ENABLE_AMR
1194 764 : std::vector<const Elem *> neighbors_offspring;
1195 : #endif
1196 :
1197 : // This is costly, and we only need to do it if the mesh has
1198 : // changed since we last partitioned... but the mesh probably has
1199 : // changed since we last partitioned, and if it hasn't we don't
1200 : // have a reliable way to be sure of that.
1201 312381 : _find_global_index_by_pid_map(mesh);
1202 :
1203 382 : dof_id_type first_local_elem = 0;
1204 1834508 : for (auto pid : make_range(mesh.processor_id()))
1205 1522318 : first_local_elem += _n_active_elem_on_proc[pid];
1206 :
1207 311999 : _dual_graph.clear();
1208 312381 : _dual_graph.resize(n_active_local_elem);
1209 312381 : _local_id_to_elem.resize(n_active_local_elem);
1210 :
1211 : // We may need to communicate constraint-row-based connections
1212 : // between processors
1213 : std::unordered_map<processor_id_type,
1214 764 : std::vector<std::pair<dof_id_type, dof_id_type>>> connections_to_push;
1215 :
1216 6205662 : for (const auto & elem : mesh.active_local_element_ptr_range())
1217 : {
1218 11826 : libmesh_assert (_global_index_by_pid_map.count(elem->id()));
1219 : const dof_id_type global_index_by_pid =
1220 2802467 : _global_index_by_pid_map[elem->id()];
1221 :
1222 2802467 : const dof_id_type local_index =
1223 11826 : global_index_by_pid - first_local_elem;
1224 11826 : libmesh_assert_less (local_index, n_active_local_elem);
1225 :
1226 23652 : std::vector<dof_id_type> & graph_row = _dual_graph[local_index];
1227 :
1228 : // Save this off to make it easy to index later
1229 2802467 : _local_id_to_elem[local_index] = const_cast<Elem*>(elem);
1230 :
1231 : // Loop over the element's neighbors. An element
1232 : // adjacency corresponds to a face neighbor
1233 14777580 : for (auto neighbor : elem->neighbor_ptr_range())
1234 : {
1235 11963287 : if (neighbor != nullptr)
1236 : {
1237 : // If the neighbor is active treat it
1238 : // as a connection
1239 50897 : if (neighbor->active())
1240 : {
1241 50770 : libmesh_assert(_global_index_by_pid_map.count(neighbor->id()));
1242 : const dof_id_type neighbor_global_index_by_pid =
1243 10833013 : _global_index_by_pid_map[neighbor->id()];
1244 :
1245 10833013 : graph_row.push_back(neighbor_global_index_by_pid);
1246 : }
1247 :
1248 : #ifdef LIBMESH_ENABLE_AMR
1249 :
1250 : // Otherwise we need to find all of the
1251 : // neighbor's children that are connected to
1252 : // us and add them
1253 : else
1254 : {
1255 : // The side of the neighbor to which
1256 : // we are connected
1257 : const unsigned int ns =
1258 44232 : neighbor->which_neighbor_am_i (elem);
1259 127 : libmesh_assert_less (ns, neighbor->n_neighbors());
1260 :
1261 : // Get all the active children (& grandchildren, etc...)
1262 : // of the neighbor
1263 :
1264 : // FIXME - this is the wrong thing, since we
1265 : // should be getting the active family tree on
1266 : // our side only. But adding too many graph
1267 : // links may cause hanging nodes to tend to be
1268 : // on partition interiors, which would reduce
1269 : // communication overhead for constraint
1270 : // equations, so we'll leave it.
1271 :
1272 44232 : neighbor->active_family_tree (neighbors_offspring);
1273 :
1274 : // Get all the neighbor's children that
1275 : // live on that side and are thus connected
1276 : // to us
1277 293509 : for (const auto & child : neighbors_offspring)
1278 : {
1279 : // This does not assume a level-1 mesh.
1280 : // Note that since children have sides numbered
1281 : // coincident with the parent then this is a sufficient test.
1282 249866 : if (child->neighbor_ptr(ns) == elem)
1283 : {
1284 254 : libmesh_assert (child->active());
1285 254 : libmesh_assert (_global_index_by_pid_map.count(child->id()));
1286 : const dof_id_type child_global_index_by_pid =
1287 118113 : _global_index_by_pid_map[child->id()];
1288 :
1289 118113 : graph_row.push_back(child_global_index_by_pid);
1290 : }
1291 : }
1292 : }
1293 :
1294 : #endif /* ifdef LIBMESH_ENABLE_AMR */
1295 :
1296 :
1297 : }
1298 : }
1299 :
1300 3695315 : if ((elem->dim() < LIBMESH_DIM) &&
1301 892848 : elem->interior_parent())
1302 : {
1303 : // get all relevant interior elements
1304 0 : std::set<const Elem *> neighbor_set;
1305 896 : elem->find_interior_neighbors(neighbor_set);
1306 :
1307 1792 : for (const auto & neighbor : neighbor_set)
1308 : {
1309 : const dof_id_type neighbor_global_index_by_pid =
1310 896 : _global_index_by_pid_map[neighbor->id()];
1311 :
1312 896 : graph_row.push_back(neighbor_global_index_by_pid);
1313 : }
1314 : }
1315 :
1316 : // Check for any boundary neighbors
1317 2803363 : for (const auto & pr : as_range(interior_to_boundary_map.equal_range(elem)))
1318 : {
1319 896 : const Elem * neighbor = pr.second;
1320 :
1321 : const dof_id_type neighbor_global_index_by_pid =
1322 896 : _global_index_by_pid_map[neighbor->id()];
1323 :
1324 896 : graph_row.push_back(neighbor_global_index_by_pid);
1325 : }
1326 :
1327 : // Check for any constraining elements
1328 2802467 : if (!mesh_constrained_nodes.empty()) // quick test for non-IGA cases
1329 : {
1330 0 : const auto end_it = mesh_constrained_nodes.end();
1331 :
1332 : // Use a set to avoid duplicates. Use a well-defined method
1333 : // of ordering that set to make debugging easier.
1334 0 : std::set<const Elem *, CompareElemIdsByLevel> constraining_elems;
1335 146202 : for (const Node & node : elem->node_ref_range())
1336 : {
1337 118363 : if (const auto row_it = mesh_constrained_nodes.find(&node);
1338 0 : row_it != end_it)
1339 555100 : for (const auto & [pr, coef] : row_it->second)
1340 : {
1341 0 : libmesh_ignore(coef); // avoid gcc 7 warning
1342 451304 : constraining_elems.insert(pr.first);
1343 : }
1344 : }
1345 240695 : for (const Elem * constraining_elem : constraining_elems)
1346 : {
1347 : const dof_id_type constraining_global_index_by_pid =
1348 212856 : _global_index_by_pid_map[constraining_elem->id()];
1349 :
1350 212856 : graph_row.push_back(constraining_global_index_by_pid);
1351 :
1352 : // We can't be sure if the constraining element's owner sees
1353 : // the assembly element, so to get a symmetric connectivity
1354 : // graph we'll need to tell them about us to be safe.
1355 212856 : if (constraining_elem->processor_id() != mesh.processor_id())
1356 237014 : connections_to_push[constraining_elem->processor_id()].emplace_back
1357 118507 : (global_index_by_pid, constraining_global_index_by_pid);
1358 : }
1359 : }
1360 :
1361 : // Check for any constrained elements
1362 2978840 : for (const auto & pr : as_range(elems_constrained_by.equal_range(elem)))
1363 : {
1364 176373 : const Elem * constrained = pr.second;
1365 : const dof_id_type constrained_global_index_by_pid =
1366 176373 : _global_index_by_pid_map[constrained->id()];
1367 :
1368 176373 : graph_row.push_back(constrained_global_index_by_pid);
1369 :
1370 : // We can't be sure if the constrained element's owner sees
1371 : // the assembly element, so to get a symmetric connectivity
1372 : // graph we'll need to tell them about us to be safe.
1373 176373 : if (constrained->processor_id() != mesh.processor_id())
1374 164048 : connections_to_push[constrained->processor_id()].emplace_back
1375 82024 : (global_index_by_pid, constrained_global_index_by_pid);
1376 : }
1377 311617 : }
1378 :
1379 : // Partitioners like Parmetis require a symmetric adjacency matrix,
1380 : // but if we have an assembly element constrained by a spline
1381 : // NodeElem owned by another processor, it's possible that that
1382 : // processor doesn't see our assembly element. Let's push those
1383 : // entries, to ensure they're counted on both sides.
1384 : auto symmetrize_entries =
1385 4059 : [this, first_local_elem]
1386 : (processor_id_type /*src_pid*/,
1387 200531 : const std::vector<std::pair<dof_id_type, dof_id_type>> & incoming_entries)
1388 : {
1389 204590 : for (auto [i, j] : incoming_entries)
1390 : {
1391 0 : libmesh_assert_greater_equal(j, first_local_elem);
1392 200531 : const std::size_t jl = j - first_local_elem;
1393 0 : libmesh_assert_less(jl, _dual_graph.size());
1394 0 : std::vector<dof_id_type> & graph_row = _dual_graph[jl];
1395 200531 : if (std::find(graph_row.begin(), graph_row.end(), i) == graph_row.end())
1396 : {
1397 : // std::cerr << "Pushing back (" << j << ", " << i << ") from " << src_pid << std::endl;
1398 36483 : graph_row.push_back(i);
1399 : }
1400 : }
1401 311999 : };
1402 :
1403 312381 : Parallel::push_parallel_vector_data(mesh.comm(),
1404 : connections_to_push,
1405 : symmetrize_entries);
1406 :
1407 : // *Now* we should have a symmetric adjacency matrix. Let's check
1408 : // that, so any failures get caught before they e.g. confuse
1409 : // Parmetis in hard-to-debug ways. That's a global communication so
1410 : // we'll only do it in debug mode.
1411 : #ifdef DEBUG
1412 382 : auto n_proc = mesh.n_processors();
1413 764 : std::vector<dof_id_type> first_local_index_on_proc(n_proc, 0);
1414 764 : for (auto pid : make_range(1u,n_proc))
1415 382 : first_local_index_on_proc[pid] = first_local_index_on_proc[pid-1] + _n_active_elem_on_proc[pid-1];
1416 :
1417 382 : libmesh_assert_equal_to(first_local_index_on_proc[mesh.processor_id()],
1418 : first_local_elem);
1419 :
1420 764 : std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, dof_id_type>>> entries_to_send;
1421 12208 : for (auto il : index_range(_dual_graph))
1422 : {
1423 11826 : const std::vector<dof_id_type> & graph_row = _dual_graph[il];
1424 :
1425 11826 : const auto i = il + first_local_elem;
1426 :
1427 62850 : for (auto j : graph_row)
1428 : {
1429 : // Stupid graph rows aren't sorted yet...
1430 51024 : processor_id_type target_pid = 0;
1431 127340 : while (target_pid+1 < n_proc &&
1432 51024 : j >= first_local_index_on_proc[target_pid+1])
1433 25292 : ++target_pid;
1434 :
1435 51024 : entries_to_send[target_pid].emplace_back(i,j);
1436 : }
1437 : }
1438 :
1439 764 : std::vector<std::tuple<processor_id_type, dof_id_type, dof_id_type>> bad_entries;
1440 :
1441 : auto check_incoming_entries =
1442 : [this, first_local_elem, &bad_entries]
1443 : (processor_id_type src_pid,
1444 204824 : const std::vector<std::pair<dof_id_type, dof_id_type>> & incoming_entries)
1445 : {
1446 51752 : for (auto [i, j] : incoming_entries)
1447 : {
1448 51024 : if (j < first_local_elem)
1449 : {
1450 0 : bad_entries.emplace_back(src_pid,i,j);
1451 0 : continue;
1452 : }
1453 102048 : const std::size_t jl = j - first_local_elem;
1454 51024 : if (jl >= _dual_graph.size())
1455 : {
1456 0 : bad_entries.emplace_back(src_pid,i,j);
1457 0 : continue;
1458 : }
1459 51024 : const std::vector<dof_id_type> & graph_row = _dual_graph[jl];
1460 51024 : if (std::find(graph_row.begin(), graph_row.end(), i) ==
1461 102048 : graph_row.end())
1462 : {
1463 0 : bad_entries.emplace_back(src_pid,i,j);
1464 0 : continue;
1465 : }
1466 : }
1467 728 : };
1468 :
1469 : // Keep any failures in sync in parallel, for easier debugging in
1470 : // unit tests where the failure exception will be caught.
1471 : Parallel::push_parallel_vector_data
1472 382 : (mesh.comm(), entries_to_send, check_incoming_entries);
1473 382 : bool bad_entries_exist = !bad_entries.empty();
1474 382 : mesh.comm().max(bad_entries_exist);
1475 382 : if (bad_entries_exist)
1476 : {
1477 : #if 0 // Optional verbosity for if this breaks again...
1478 : if (!bad_entries.empty())
1479 : {
1480 : std::cerr << "Bad entries on processor " << mesh.processor_id() << ": ";
1481 : for (auto [p, i, j] : bad_entries)
1482 : std::cerr << '(' << p << ", " << i << ", " << j << "), ";
1483 : std::cerr << std::endl;
1484 : }
1485 : #endif
1486 0 : libmesh_error_msg("Asymmetric partitioner graph detected");
1487 : }
1488 : #endif
1489 312381 : }
1490 :
1491 54981 : void Partitioner::assign_partitioning (MeshBase & mesh, const std::vector<dof_id_type> & parts)
1492 : {
1493 882 : LOG_SCOPE("assign_partitioning()", "Partitioner");
1494 :
1495 : // This function must be run on all processors at once
1496 294 : libmesh_parallel_only(mesh.comm());
1497 :
1498 294 : dof_id_type first_local_elem = 0;
1499 277223 : for (auto pid : make_range(mesh.processor_id()))
1500 222389 : first_local_elem += _n_active_elem_on_proc[pid];
1501 :
1502 : #ifndef NDEBUG
1503 294 : const dof_id_type n_active_local_elem = mesh.n_active_local_elem();
1504 : #endif
1505 :
1506 : std::map<processor_id_type, std::vector<dof_id_type>>
1507 588 : requested_ids;
1508 :
1509 : // Results to gather from each processor - kept in a map so we
1510 : // do only one loop over elements after all receives are done.
1511 : std::map<processor_id_type, std::vector<processor_id_type>>
1512 588 : filled_request;
1513 :
1514 14363317 : for (auto & elem : mesh.active_element_ptr_range())
1515 : {
1516 : // we need to get the index from the owning processor
1517 : // (note we cannot assign it now -- we are iterating
1518 : // over elements again and this will be bad!)
1519 14272841 : requested_ids[elem->processor_id()].push_back(elem->id());
1520 54393 : }
1521 :
1522 : auto gather_functor =
1523 445731 : [this,
1524 : & parts,
1525 : #ifndef NDEBUG
1526 : & mesh,
1527 : n_active_local_elem,
1528 : #endif
1529 : first_local_elem]
1530 : (processor_id_type, const std::vector<dof_id_type> & ids,
1531 28604434 : std::vector<processor_id_type> & data)
1532 : {
1533 1176 : const std::size_t ids_size = ids.size();
1534 446907 : data.resize(ids.size());
1535 :
1536 14700556 : for (std::size_t i=0; i != ids_size; i++)
1537 : {
1538 14272841 : const dof_id_type requested_elem_index = ids[i];
1539 :
1540 19192 : libmesh_assert(_global_index_by_pid_map.count(requested_elem_index));
1541 :
1542 : const dof_id_type global_index_by_pid =
1543 14253649 : _global_index_by_pid_map[requested_elem_index];
1544 :
1545 14253649 : const dof_id_type local_index =
1546 19192 : global_index_by_pid - first_local_elem;
1547 :
1548 19192 : libmesh_assert_less (local_index, parts.size());
1549 19192 : libmesh_assert_less (local_index, n_active_local_elem);
1550 :
1551 : const processor_id_type elem_procid =
1552 14253649 : cast_int<processor_id_type>(parts[local_index]);
1553 :
1554 19192 : libmesh_assert_less (elem_procid, mesh.n_partitions());
1555 :
1556 14272841 : data[i] = elem_procid;
1557 : }
1558 501594 : };
1559 :
1560 : auto action_functor =
1561 445731 : [&filled_request]
1562 : (processor_id_type pid,
1563 : const std::vector<dof_id_type> &,
1564 446319 : const std::vector<processor_id_type> & new_procids)
1565 : {
1566 446907 : filled_request[pid] = new_procids;
1567 55569 : };
1568 :
1569 : // Trade requests with other processors
1570 294 : const processor_id_type * ex = nullptr;
1571 : Parallel::pull_parallel_vector_data
1572 54981 : (mesh.comm(), requested_ids, gather_functor, action_functor, ex);
1573 :
1574 : // and finally assign the partitioning.
1575 : // note we are iterating in exactly the same order
1576 : // used to build up the request, so we can expect the
1577 : // required entries to be in the proper sequence.
1578 55569 : std::vector<unsigned int> counters(mesh.n_processors(), 0);
1579 14363317 : for (auto & elem : mesh.active_element_ptr_range())
1580 : {
1581 14253649 : const processor_id_type current_pid = elem->processor_id();
1582 :
1583 19192 : libmesh_assert_less (counters[current_pid], requested_ids[current_pid].size());
1584 :
1585 : const processor_id_type elem_procid =
1586 14253649 : filled_request[current_pid][counters[current_pid]++];
1587 :
1588 19192 : libmesh_assert_less (elem_procid, mesh.n_partitions());
1589 14253649 : elem->processor_id() = elem_procid;
1590 54393 : }
1591 54981 : }
1592 :
1593 :
1594 : } // namespace libMesh
|
