Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
petsc_communication_pattern.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2023 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16
17#ifdef DEAL_II_WITH_PETSC
18
19# include <deal.II/base/mpi.h>
20
22
24// Shorthand notation for PETSc error codes.
25# define AssertPETSc(code) \
26 do \
27 { \
28 PetscErrorCode ierr = (code); \
29 AssertThrow(ierr == 0, ExcPETScError(ierr)); \
30 } \
31 while (false)
32
33namespace PETScWrappers
34{
38
39
40
45
46
47
48 void
50 const IndexSet &ghost_indices,
51 const MPI_Comm communicator)
52 {
53 clear();
54
55 PetscLayout layout;
56 AssertPETSc(PetscLayoutCreate(communicator, &layout));
57 AssertPETSc(PetscLayoutSetLocalSize(layout, local_size));
58 AssertPETSc(PetscLayoutSetUp(layout));
59
60 PetscInt start, end;
61 AssertPETSc(PetscLayoutGetRange(layout, &start, &end));
62
63 IndexSet want;
64 want.add_range(start, end);
65 want.add_indices(ghost_indices);
66 want.compress();
67
68 const PetscInt *idxs;
69 PetscInt n;
70 IS is = want.make_petsc_is(communicator);
71 AssertPETSc(ISGetLocalSize(is, &n));
72 AssertPETSc(ISGetIndices(is, &idxs));
73
74 AssertPETSc(PetscSFCreate(communicator, &sf));
76 PetscSFSetGraphLayout(sf, layout, n, nullptr, PETSC_OWN_POINTER, idxs));
77 AssertPETSc(PetscSFSetUp(sf));
78
79 AssertPETSc(ISRestoreIndices(is, &idxs));
80 AssertPETSc(ISDestroy(&is));
81 AssertPETSc(PetscLayoutDestroy(&layout));
82 }
83
84
85
86 void
87 CommunicationPattern::reinit(const IndexSet &locally_owned_indices,
88 const IndexSet &ghost_indices,
89 const MPI_Comm communicator)
90 {
91 std::vector<types::global_dof_index> in_deal;
92 locally_owned_indices.fill_index_vector(in_deal);
93 std::vector<PetscInt> in_petsc(in_deal.begin(), in_deal.end());
94
95 std::vector<types::global_dof_index> out_deal;
96 ghost_indices.fill_index_vector(out_deal);
97 std::vector<PetscInt> out_petsc(out_deal.begin(), out_deal.end());
98
99 std::vector<PetscInt> dummy;
100
101 this->do_reinit(in_petsc, dummy, out_petsc, dummy, communicator);
102 }
103
104
105
106 void
108 const std::vector<types::global_dof_index> &indices_has,
109 const std::vector<types::global_dof_index> &indices_want,
110 const MPI_Comm communicator)
111 {
112 // Clean vectors from numbers::invalid_dof_index (indicating padding)
113 std::vector<PetscInt> indices_has_clean, indices_has_loc;
114 std::vector<PetscInt> indices_want_clean, indices_want_loc;
115 indices_want_clean.reserve(indices_want.size());
116 indices_want_loc.reserve(indices_want.size());
117 indices_has_clean.reserve(indices_has.size());
118 indices_has_loc.reserve(indices_has.size());
119
120 PetscInt loc = 0;
121 bool has_invalid = false;
122 for (const auto i : indices_has)
123 {
125 {
126 indices_has_clean.push_back(static_cast<PetscInt>(i));
127 indices_has_loc.push_back(loc);
128 }
129 else
130 has_invalid = true;
131 ++loc;
132 }
133 if (!has_invalid)
134 indices_has_loc.clear();
135
136 loc = 0;
137 has_invalid = false;
138 for (const auto i : indices_want)
139 {
141 {
142 indices_want_clean.push_back(static_cast<PetscInt>(i));
143 indices_want_loc.push_back(loc);
144 }
145 else
146 has_invalid = true;
147 ++loc;
148 }
149 if (!has_invalid)
150 indices_want_loc.clear();
151
152 this->do_reinit(indices_has_clean,
153 indices_has_loc,
154 indices_want_clean,
155 indices_want_loc,
156 communicator);
157 }
158
159
160
161 void
162 CommunicationPattern::do_reinit(const std::vector<PetscInt> &inidx,
163 const std::vector<PetscInt> &inloc,
164 const std::vector<PetscInt> &outidx,
165 const std::vector<PetscInt> &outloc,
166 const MPI_Comm communicator)
167 {
168 clear();
169
170 // inidx is assumed to be unstructured and non-overlapping.
171 // However, it may have holes in it and not be a full cover.
172 //
173 // We create two PETSc SFs and compose them to get
174 // the final communication pattern
175 //
176 // sf1 : local distributed to tmp
177 // sf2 : tmp to local with ghosts
178 // sf(x) = sf2(sf1(x))
179 PetscSF sf1, sf2;
180
181 // First create an SF where leaves are inidx (at location inloc)
182 // and roots are unique indices in contiguous way
183 // Code adapted from MatZeroRowsMapLocal_Private in PETSc
184 PetscInt n = static_cast<PetscInt>(inidx.size());
185 PetscInt lN = n > 0 ? *std::max_element(inidx.begin(), inidx.end()) : -1;
186 PetscInt N, nl;
187
188 Utilities::MPI::internal::all_reduce<PetscInt>(
189 MPI_MAX,
191 communicator,
192 ArrayView<PetscInt>(&N, 1));
193
194 PetscSFNode *remotes;
195 AssertPETSc(PetscMalloc1(n, &remotes));
196
197 PetscLayout layout;
198 AssertPETSc(PetscLayoutCreate(communicator, &layout));
199 AssertPETSc(PetscLayoutSetSize(layout, N + 1));
200 AssertPETSc(PetscLayoutSetUp(layout));
201 AssertPETSc(PetscLayoutGetLocalSize(layout, &nl));
202
203 const PetscInt *ranges;
204 AssertPETSc(PetscLayoutGetRanges(layout, &ranges));
205
206 PetscInt cnt = 0;
207# if DEAL_II_PETSC_VERSION_GTE(3, 13, 0)
208 PetscMPIInt owner = 0;
209# else
210 PetscInt owner = 0;
211# endif
212 for (const auto idx : inidx)
213 {
214 // short-circuit the search if the last owner owns this index too
215 if (idx < ranges[owner] || ranges[owner + 1] <= idx)
216 {
217 AssertPETSc(PetscLayoutFindOwner(layout, idx, &owner));
218 }
219 remotes[cnt].rank = owner;
220 remotes[cnt].index = idx - ranges[owner];
221 ++cnt;
222 }
223
224 AssertPETSc(PetscSFCreate(communicator, &sf2));
225 AssertPETSc(PetscSFSetGraph(sf2,
226 nl,
227 n,
228 const_cast<PetscInt *>(
229 inloc.size() > 0 ? inloc.data() : nullptr),
230 PETSC_COPY_VALUES,
231 remotes,
232 PETSC_OWN_POINTER));
233 AssertPETSc(PetscSFSetUp(sf2));
234 // We need to invert root and leaf space to create the first SF
235 AssertPETSc(PetscSFCreateInverseSF(sf2, &sf1));
236 AssertPETSc(PetscSFDestroy(&sf2));
237
238 // Now create the SF from the contiguous space to the local output space
239 n = static_cast<PetscInt>(outidx.size());
240 AssertPETSc(PetscSFCreate(communicator, &sf2));
241 AssertPETSc(PetscSFSetGraphLayout(
242 sf2,
243 layout,
244 n,
245 const_cast<PetscInt *>(outloc.size() > 0 ? outloc.data() : nullptr),
246 PETSC_COPY_VALUES,
247 const_cast<PetscInt *>(n > 0 ? outidx.data() : nullptr)));
248 AssertPETSc(PetscSFSetUp(sf2));
249
250 // The final SF is the composition of the two
251 AssertPETSc(PetscSFCompose(sf1, sf2, &sf));
252
253 // Cleanup
254 AssertPETSc(PetscLayoutDestroy(&layout));
255 AssertPETSc(PetscSFDestroy(&sf1));
256 AssertPETSc(PetscSFDestroy(&sf2));
257 }
258
259
260
261 void
263 {
264 AssertPETSc(PetscSFDestroy(&sf));
265 }
266
267
268
271 {
272 return PetscObjectComm(reinterpret_cast<PetscObject>(sf));
273 }
274
275
276
277 template <typename Number>
278 void
280 const ArrayView<const Number> &src,
281 const ArrayView<Number> &dst) const
282 {
284
285# if DEAL_II_PETSC_VERSION_LT(3, 15, 0)
286 AssertPETSc(PetscSFBcastBegin(sf, datatype, src.data(), dst.data()));
287# else
289 PetscSFBcastBegin(sf, datatype, src.data(), dst.data(), MPI_REPLACE));
290# endif
291 }
292
293
294
295 template <typename Number>
296 void
298 const ArrayView<const Number> &src,
299 const ArrayView<Number> &dst) const
300 {
302
303# if DEAL_II_PETSC_VERSION_LT(3, 15, 0)
304 AssertPETSc(PetscSFBcastEnd(sf, datatype, src.data(), dst.data()));
305# else
307 PetscSFBcastEnd(sf, datatype, src.data(), dst.data(), MPI_REPLACE));
308# endif
309 }
310
311
312
313 template <typename Number>
314 void
322
323
324
325 template <typename Number>
326 void
329 const ArrayView<const Number> &src,
330 const ArrayView<Number> &dst) const
331 {
332 MPI_Op mpiop = (op == VectorOperation::insert) ? MPI_REPLACE : MPI_SUM;
334
336 PetscSFReduceBegin(sf, datatype, src.data(), dst.data(), mpiop));
337 }
338
339
340
341 template <typename Number>
342 void
345 const ArrayView<const Number> &src,
346 const ArrayView<Number> &dst) const
347 {
348 MPI_Op mpiop = (op == VectorOperation::insert) ? MPI_REPLACE : MPI_SUM;
350
351 AssertPETSc(PetscSFReduceEnd(sf, datatype, src.data(), dst.data(), mpiop));
352 }
353
354
355
356 template <typename Number>
357 void
366
367# ifndef DOXYGEN
368
369 // Partitioner
370
372 : ghost()
373 , larger_ghost()
374 , ghost_indices_data()
375 , n_ghost_indices_data(numbers::invalid_dof_index)
376 , n_ghost_indices_larger(numbers::invalid_dof_index)
377 {}
378
379
380
381 void
382 Partitioner::reinit(const IndexSet &locally_owned_indices,
383 const IndexSet &ghost_indices,
384 const MPI_Comm communicator)
385 {
387 ghost_indices_data.subtract_set(locally_owned_indices);
389
390 ghost.reinit(locally_owned_indices, ghost_indices_data, communicator);
392
395 }
396
397 void
398 Partitioner::reinit(const IndexSet &locally_owned_indices,
399 const IndexSet &ghost_indices,
400 const IndexSet &larger_ghost_indices,
401 const MPI_Comm communicator)
402 {
403 std::vector<types::global_dof_index> local_indices;
404 locally_owned_indices.fill_index_vector(local_indices);
405
407 ghost_indices_data.subtract_set(locally_owned_indices);
409
410 std::vector<types::global_dof_index> expanded_ghost_indices(
411 larger_ghost_indices.n_elements(), numbers::invalid_dof_index);
412 for (auto index : ghost_indices_data)
413 {
414 Assert(larger_ghost_indices.is_element(index),
415 ExcMessage("The given larger ghost index set must contain "
416 "all indices in the actual index set."));
417 auto tmp_index = larger_ghost_indices.index_within_set(index);
418 expanded_ghost_indices[tmp_index] = index;
419 }
420
421 ghost.reinit(locally_owned_indices, ghost_indices_data, communicator);
422 larger_ghost.reinit(local_indices, expanded_ghost_indices, communicator);
424 n_ghost_indices_larger = larger_ghost_indices.n_elements();
425 }
426
429 {
431 }
432
433 template <typename Number>
434 void
436 const ArrayView<Number> &dst) const
437 {
438 if (dst.size() == n_ghost_indices_larger)
439 {
441 }
442 else
443 {
445 }
446 }
447
448 template <typename Number>
449 void
451 const ArrayView<const Number> &src,
452 const ArrayView<Number> &dst) const
453 {
454 if (dst.size() == n_ghost_indices_larger)
455 {
457 }
458 else
459 {
461 }
462 }
463
464 template <typename Number>
465 void
467 const ArrayView<Number> &dst) const
468 {
471 }
472
473 template <typename Number>
474 void
477 const ArrayView<const Number> &src,
478 const ArrayView<Number> &dst) const
479 {
480 if (src.size() == n_ghost_indices_larger)
481 {
483 }
484 else
485 {
487 }
488 }
489
490 template <typename Number>
491 void
494 const ArrayView<const Number> &src,
495 const ArrayView<Number> &dst) const
496 {
497 if (src.size() == n_ghost_indices_larger)
498 {
500 }
501 else
502 {
504 }
505 }
506
507 template <typename Number>
508 void
510 const ArrayView<const Number> &src,
511 const ArrayView<Number> &dst) const
512 {
515 }
516# endif
517} // namespace PETScWrappers
518
519// Explicit instantiations
520# include "petsc_communication_pattern.inst"
521
523
524#endif // DEAL_II_WITH_PETSC
value_type * data() const noexcept
Definition array_view.h:661
std::size_t size() const
Definition array_view.h:684
IS make_petsc_is(const MPI_Comm communicator=MPI_COMM_WORLD) const
void fill_index_vector(std::vector< size_type > &indices) const
Definition index_set.cc:945
size_type index_within_set(const size_type global_index) const
Definition index_set.h:2001
size_type n_elements() const
Definition index_set.h:1934
bool is_element(const size_type index) const
Definition index_set.h:1894
void subtract_set(const IndexSet &other)
Definition index_set.cc:473
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1803
void compress() const
Definition index_set.h:1784
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1831
void import_from_ghosted_array(const VectorOperation::values op, const ArrayView< const Number > &ghost_array, const ArrayView< Number > &locally_owned_array) const
void export_to_ghosted_array_start(const ArrayView< const Number > &locally_owned_array, const ArrayView< Number > &ghost_array) const
void export_to_ghosted_array(const ArrayView< const Number > &locally_owned_array, const ArrayView< Number > &ghost_array) const
void export_to_ghosted_array_finish(const ArrayView< const Number > &locally_owned_array, const ArrayView< Number > &ghost_array) const
void import_from_ghosted_array_finish(const VectorOperation::values op, const ArrayView< const Number > &ghost_array, const ArrayView< Number > &locally_owned_array) const
void do_reinit(const std::vector< PetscInt > &inidx, const std::vector< PetscInt > &inloc, const std::vector< PetscInt > &outidx, const std::vector< PetscInt > &outloc, const MPI_Comm communicator)
virtual void reinit(const IndexSet &locally_owned_indices, const IndexSet &ghost_indices, const MPI_Comm communicator) override
void import_from_ghosted_array_start(const VectorOperation::values op, const ArrayView< const Number > &ghost_array, const ArrayView< Number > &locally_owned_array) const
types::global_dof_index n_ghost_indices_data
virtual void reinit(const IndexSet &locally_owned_indices, const IndexSet &ghost_indices, const MPI_Comm communicator) override
void import_from_ghosted_array(const VectorOperation::values op, const ArrayView< const Number > &ghost_array, const ArrayView< Number > &locally_owned_array) const
void export_to_ghosted_array_finish(const ArrayView< const Number > &locally_owned_array, const ArrayView< Number > &ghost_array) const
MPI_Comm get_mpi_communicator() const override
void export_to_ghosted_array(const ArrayView< const Number > &locally_owned_array, const ArrayView< Number > &ghost_array) const
void export_to_ghosted_array_start(const ArrayView< const Number > &locally_owned_array, const ArrayView< Number > &ghost_array) const
void import_from_ghosted_array_start(const VectorOperation::values op, const ArrayView< const Number > &ghost_array, const ArrayView< Number > &locally_owned_array) const
void import_from_ghosted_array_finish(const VectorOperation::values op, const ArrayView< const Number > &ghost_array, const ArrayView< Number > &locally_owned_array) const
types::global_dof_index n_ghost_indices_larger
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
const MPI_Datatype mpi_type_id_for_type
Definition mpi.h:1674
const types::global_dof_index invalid_dof_index
Definition types.h:252
#define AssertPETSc(code)