deal.II version GIT relicensing-2017-g7ae82fad8b 2024-10-22 19:20:00+00:00
\(\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 AssertIndexRange(local_size, ghost_indices.size() + 1);
55 // If the size of the index set can be converted to a PetscInt then every
56 // index can also be converted
57 AssertThrowIntegerConversion(ghost_indices.size(),
58 static_cast<PetscInt>(ghost_indices.size()));
59
60 PetscLayout layout;
61 AssertPETSc(PetscLayoutCreate(communicator, &layout));
62 const auto petsc_local_size = static_cast<PetscInt>(local_size);
63 AssertThrowIntegerConversion(local_size, petsc_local_size);
64 AssertPETSc(PetscLayoutSetLocalSize(layout, petsc_local_size));
65 AssertPETSc(PetscLayoutSetUp(layout));
66
67 PetscInt start, end;
68 AssertPETSc(PetscLayoutGetRange(layout, &start, &end));
69
70 IndexSet want;
71 want.add_range(start, end);
72 want.add_indices(ghost_indices);
73 want.compress();
74
75 const PetscInt *idxs;
76 PetscInt n;
77 IS is = want.make_petsc_is(communicator);
78 AssertPETSc(ISGetLocalSize(is, &n));
79 AssertPETSc(ISGetIndices(is, &idxs));
80
81 AssertPETSc(PetscSFCreate(communicator, &sf));
83 PetscSFSetGraphLayout(sf, layout, n, nullptr, PETSC_OWN_POINTER, idxs));
84 AssertPETSc(PetscSFSetUp(sf));
85
86 AssertPETSc(ISRestoreIndices(is, &idxs));
87 AssertPETSc(ISDestroy(&is));
88 AssertPETSc(PetscLayoutDestroy(&layout));
89 }
90
91
92
93 void
94 CommunicationPattern::reinit(const IndexSet &locally_owned_indices,
95 const IndexSet &ghost_indices,
96 const MPI_Comm communicator)
97 {
98 // If the sizes of the index sets can be converted to PetscInts then every
99 // index can also be converted
100 AssertThrowIntegerConversion(static_cast<PetscInt>(
101 locally_owned_indices.size()),
102 locally_owned_indices.size());
103 AssertThrowIntegerConversion(static_cast<PetscInt>(ghost_indices.size()),
104 ghost_indices.size());
105
106 const auto in_deal = locally_owned_indices.get_index_vector();
107 std::vector<PetscInt> in_petsc(in_deal.begin(), in_deal.end());
108
109 const auto out_deal = ghost_indices.get_index_vector();
110 std::vector<PetscInt> out_petsc(out_deal.begin(), out_deal.end());
111
112 std::vector<PetscInt> dummy;
113
114 this->do_reinit(in_petsc, dummy, out_petsc, dummy, communicator);
115 }
116
117
118
119 void
121 const std::vector<types::global_dof_index> &indices_has,
122 const std::vector<types::global_dof_index> &indices_want,
123 const MPI_Comm communicator)
124 {
125 // Clean vectors from numbers::invalid_dof_index (indicating padding)
126 std::vector<PetscInt> indices_has_clean, indices_has_loc;
127 std::vector<PetscInt> indices_want_clean, indices_want_loc;
128 indices_want_clean.reserve(indices_want.size());
129 indices_want_loc.reserve(indices_want.size());
130 indices_has_clean.reserve(indices_has.size());
131 indices_has_loc.reserve(indices_has.size());
132
133 PetscInt loc = 0;
134 bool has_invalid = false;
135 for (const auto i : indices_has)
136 {
138 {
139 const auto petsc_i = static_cast<PetscInt>(i);
141 indices_has_clean.push_back(petsc_i);
142 indices_has_loc.push_back(loc);
143 }
144 else
145 has_invalid = true;
146 ++loc;
147 }
148 if (!has_invalid)
149 indices_has_loc.clear();
150
151 loc = 0;
152 has_invalid = false;
153 for (const auto i : indices_want)
154 {
156 {
157 const auto petsc_i = static_cast<PetscInt>(i);
159 indices_want_clean.push_back(petsc_i);
160 indices_want_loc.push_back(loc);
161 }
162 else
163 has_invalid = true;
164 ++loc;
165 }
166 if (!has_invalid)
167 indices_want_loc.clear();
168
169 this->do_reinit(indices_has_clean,
170 indices_has_loc,
171 indices_want_clean,
172 indices_want_loc,
173 communicator);
174 }
175
176
177
178 void
179 CommunicationPattern::do_reinit(const std::vector<PetscInt> &inidx,
180 const std::vector<PetscInt> &inloc,
181 const std::vector<PetscInt> &outidx,
182 const std::vector<PetscInt> &outloc,
183 const MPI_Comm communicator)
184 {
185 clear();
186
187 // inidx is assumed to be unstructured and non-overlapping.
188 // However, it may have holes in it and not be a full cover.
189 //
190 // We create two PETSc SFs and compose them to get
191 // the final communication pattern
192 //
193 // sf1 : local distributed to tmp
194 // sf2 : tmp to local with ghosts
195 // sf(x) = sf2(sf1(x))
196 PetscSF sf1, sf2;
197
198 // First create an SF where leaves are inidx (at location inloc)
199 // and roots are unique indices in contiguous way
200 // Code adapted from MatZeroRowsMapLocal_Private in PETSc
201 PetscInt n = static_cast<PetscInt>(inidx.size());
202 PetscInt lN = n > 0 ? *std::max_element(inidx.begin(), inidx.end()) : -1;
203 PetscInt N, nl;
204
205 Utilities::MPI::internal::all_reduce<PetscInt>(
206 MPI_MAX,
208 communicator,
209 ArrayView<PetscInt>(&N, 1));
210
211 PetscSFNode *remotes;
212 AssertPETSc(PetscMalloc1(n, &remotes));
213
214 PetscLayout layout;
215 AssertPETSc(PetscLayoutCreate(communicator, &layout));
216 AssertPETSc(PetscLayoutSetSize(layout, N + 1));
217 AssertPETSc(PetscLayoutSetUp(layout));
218 AssertPETSc(PetscLayoutGetLocalSize(layout, &nl));
219
220 const PetscInt *ranges;
221 AssertPETSc(PetscLayoutGetRanges(layout, &ranges));
222
223 PetscInt cnt = 0;
224# if DEAL_II_PETSC_VERSION_GTE(3, 13, 0)
225 PetscMPIInt owner = 0;
226# else
227 PetscInt owner = 0;
228# endif
229 for (const auto idx : inidx)
230 {
231 // short-circuit the search if the last owner owns this index too
232 if (idx < ranges[owner] || ranges[owner + 1] <= idx)
233 {
234 AssertPETSc(PetscLayoutFindOwner(layout, idx, &owner));
235 }
236 remotes[cnt].rank = owner;
237 remotes[cnt].index = idx - ranges[owner];
238 ++cnt;
239 }
240
241 AssertPETSc(PetscSFCreate(communicator, &sf2));
242 AssertPETSc(PetscSFSetGraph(sf2,
243 nl,
244 n,
245 const_cast<PetscInt *>(
246 inloc.size() > 0 ? inloc.data() : nullptr),
247 PETSC_COPY_VALUES,
248 remotes,
249 PETSC_OWN_POINTER));
250 AssertPETSc(PetscSFSetUp(sf2));
251 // We need to invert root and leaf space to create the first SF
252 AssertPETSc(PetscSFCreateInverseSF(sf2, &sf1));
253 AssertPETSc(PetscSFDestroy(&sf2));
254
255 // Now create the SF from the contiguous space to the local output space
256 n = static_cast<PetscInt>(outidx.size());
257 AssertPETSc(PetscSFCreate(communicator, &sf2));
258 AssertPETSc(PetscSFSetGraphLayout(
259 sf2,
260 layout,
261 n,
262 const_cast<PetscInt *>(outloc.size() > 0 ? outloc.data() : nullptr),
263 PETSC_COPY_VALUES,
264 const_cast<PetscInt *>(n > 0 ? outidx.data() : nullptr)));
265 AssertPETSc(PetscSFSetUp(sf2));
266
267 // The final SF is the composition of the two
268 AssertPETSc(PetscSFCompose(sf1, sf2, &sf));
269
270 // Cleanup
271 AssertPETSc(PetscLayoutDestroy(&layout));
272 AssertPETSc(PetscSFDestroy(&sf1));
273 AssertPETSc(PetscSFDestroy(&sf2));
274 }
275
276
277
278 void
280 {
281 AssertPETSc(PetscSFDestroy(&sf));
282 }
283
284
285
288 {
289 return PetscObjectComm(reinterpret_cast<PetscObject>(sf));
290 }
291
292
293
294 template <typename Number>
295 void
297 const ArrayView<const Number> &src,
298 const ArrayView<Number> &dst) const
299 {
300 auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
301
302# if DEAL_II_PETSC_VERSION_LT(3, 15, 0)
303 AssertPETSc(PetscSFBcastBegin(sf, datatype, src.data(), dst.data()));
304# else
306 PetscSFBcastBegin(sf, datatype, src.data(), dst.data(), MPI_REPLACE));
307# endif
308 }
309
310
311
312 template <typename Number>
313 void
315 const ArrayView<const Number> &src,
316 const ArrayView<Number> &dst) const
317 {
318 auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
319
320# if DEAL_II_PETSC_VERSION_LT(3, 15, 0)
321 AssertPETSc(PetscSFBcastEnd(sf, datatype, src.data(), dst.data()));
322# else
324 PetscSFBcastEnd(sf, datatype, src.data(), dst.data(), MPI_REPLACE));
325# endif
326 }
327
328
329
330 template <typename Number>
331 void
339
340
341
342 template <typename Number>
343 void
346 const ArrayView<const Number> &src,
347 const ArrayView<Number> &dst) const
348 {
349 MPI_Op mpiop = (op == VectorOperation::insert) ? MPI_REPLACE : MPI_SUM;
350 auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
351
353 PetscSFReduceBegin(sf, datatype, src.data(), dst.data(), mpiop));
354 }
355
356
357
358 template <typename Number>
359 void
362 const ArrayView<const Number> &src,
363 const ArrayView<Number> &dst) const
364 {
365 MPI_Op mpiop = (op == VectorOperation::insert) ? MPI_REPLACE : MPI_SUM;
366 auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
367
368 AssertPETSc(PetscSFReduceEnd(sf, datatype, src.data(), dst.data(), mpiop));
369 }
370
371
372
373 template <typename Number>
374 void
383
384# ifndef DOXYGEN
385
386 // Partitioner
387
389 : ghost()
390 , larger_ghost()
391 , ghost_indices_data()
392 , n_ghost_indices_data(numbers::invalid_dof_index)
393 , n_ghost_indices_larger(numbers::invalid_dof_index)
394 {}
395
396
397
398 void
399 Partitioner::reinit(const IndexSet &locally_owned_indices,
400 const IndexSet &ghost_indices,
401 const MPI_Comm communicator)
402 {
404 ghost_indices_data.subtract_set(locally_owned_indices);
406
407 ghost.reinit(locally_owned_indices, ghost_indices_data, communicator);
409
412 }
413
414 void
415 Partitioner::reinit(const IndexSet &locally_owned_indices,
416 const IndexSet &ghost_indices,
417 const IndexSet &larger_ghost_indices,
418 const MPI_Comm communicator)
419 {
421 ghost_indices_data.subtract_set(locally_owned_indices);
423
424 std::vector<types::global_dof_index> expanded_ghost_indices(
425 larger_ghost_indices.n_elements(), numbers::invalid_dof_index);
426 for (auto index : ghost_indices_data)
427 {
428 Assert(larger_ghost_indices.is_element(index),
429 ExcMessage("The given larger ghost index set must contain "
430 "all indices in the actual index set."));
431 auto tmp_index = larger_ghost_indices.index_within_set(index);
432 expanded_ghost_indices[tmp_index] = index;
433 }
434
435 ghost.reinit(locally_owned_indices, ghost_indices_data, communicator);
436 larger_ghost.reinit(locally_owned_indices.get_index_vector(),
437 expanded_ghost_indices,
438 communicator);
440 n_ghost_indices_larger = larger_ghost_indices.n_elements();
441 }
442
445 {
447 }
448
449 template <typename Number>
450 void
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<const Number> &src,
468 const ArrayView<Number> &dst) const
469 {
470 if (dst.size() == n_ghost_indices_larger)
471 {
473 }
474 else
475 {
477 }
478 }
479
480 template <typename Number>
481 void
483 const ArrayView<Number> &dst) const
484 {
487 }
488
489 template <typename Number>
490 void
493 const ArrayView<const Number> &src,
494 const ArrayView<Number> &dst) const
495 {
496 if (src.size() == n_ghost_indices_larger)
497 {
499 }
500 else
501 {
503 }
504 }
505
506 template <typename Number>
507 void
510 const ArrayView<const Number> &src,
511 const ArrayView<Number> &dst) const
512 {
513 if (src.size() == n_ghost_indices_larger)
514 {
516 }
517 else
518 {
520 }
521 }
522
523 template <typename Number>
524 void
526 const ArrayView<const Number> &src,
527 const ArrayView<Number> &dst) const
528 {
531 }
532# endif
533} // namespace PETScWrappers
534
535// Explicit instantiations
536# include "petsc_communication_pattern.inst"
537
539
540#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
size_type size() const
Definition index_set.h:1776
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
std::vector< size_type > get_index_vector() const
Definition index_set.cc:926
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrowIntegerConversion(index1, index2)
const types::global_dof_index invalid_dof_index
Definition types.h:252
#define AssertPETSc(code)