Reference documentation for deal.II version 9.5.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// Copyright (C) 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17
18#ifdef DEAL_II_WITH_PETSC
19
20# include <deal.II/base/mpi.h>
21
23
25// Shorthand notation for PETSc error codes.
26# define AssertPETSc(code) \
27 do \
28 { \
29 PetscErrorCode ierr = (code); \
30 AssertThrow(ierr == 0, ExcPETScError(ierr)); \
31 } \
32 while (0)
33
34namespace PETScWrappers
35{
37 : sf(nullptr)
38 {}
39
41 {
42 clear();
43 }
44
45 void
47 const IndexSet & ghost_indices,
48 const MPI_Comm communicator)
49 {
50 clear();
51
52 PetscLayout layout;
53 AssertPETSc(PetscLayoutCreate(communicator, &layout));
54 AssertPETSc(PetscLayoutSetLocalSize(layout, local_size));
55 AssertPETSc(PetscLayoutSetUp(layout));
56
57 PetscInt start, end;
58 AssertPETSc(PetscLayoutGetRange(layout, &start, &end));
59
60 IndexSet want;
61 want.add_range(start, end);
62 want.add_indices(ghost_indices);
63 want.compress();
64
65 const PetscInt *idxs;
66 PetscInt n;
67 IS is = want.make_petsc_is(communicator);
68 AssertPETSc(ISGetLocalSize(is, &n));
69 AssertPETSc(ISGetIndices(is, &idxs));
70
71 AssertPETSc(PetscSFCreate(communicator, &sf));
73 PetscSFSetGraphLayout(sf, layout, n, nullptr, PETSC_OWN_POINTER, idxs));
74 AssertPETSc(PetscSFSetUp(sf));
75
76 AssertPETSc(ISRestoreIndices(is, &idxs));
77 AssertPETSc(ISDestroy(&is));
78 AssertPETSc(PetscLayoutDestroy(&layout));
79 }
80
81 void
82 CommunicationPattern::reinit(const IndexSet &locally_owned_indices,
83 const IndexSet &ghost_indices,
84 const MPI_Comm communicator)
85 {
86 std::vector<types::global_dof_index> in_deal;
87 locally_owned_indices.fill_index_vector(in_deal);
88 std::vector<PetscInt> in_petsc(in_deal.begin(), in_deal.end());
89
90 std::vector<types::global_dof_index> out_deal;
91 ghost_indices.fill_index_vector(out_deal);
92 std::vector<PetscInt> out_petsc(out_deal.begin(), out_deal.end());
93
94 std::vector<PetscInt> dummy;
95
96 this->do_reinit(in_petsc, dummy, out_petsc, dummy, communicator);
97 }
98
99 void
101 const std::vector<types::global_dof_index> &indices_has,
102 const std::vector<types::global_dof_index> &indices_want,
103 const MPI_Comm communicator)
104 {
105 // Clean vectors from numbers::invalid_dof_index (indicating padding)
106 std::vector<PetscInt> indices_has_clean, indices_has_loc;
107 std::vector<PetscInt> indices_want_clean, indices_want_loc;
108 indices_want_clean.reserve(indices_want.size());
109 indices_want_loc.reserve(indices_want.size());
110 indices_has_clean.reserve(indices_has.size());
111 indices_has_loc.reserve(indices_has.size());
112
113 PetscInt loc = 0;
114 bool has_invalid = false;
115 for (const auto i : indices_has)
116 {
118 {
119 indices_has_clean.push_back(static_cast<PetscInt>(i));
120 indices_has_loc.push_back(loc);
121 }
122 else
123 has_invalid = true;
124 loc++;
125 }
126 if (!has_invalid)
127 indices_has_loc.clear();
128
129 loc = 0;
130 has_invalid = false;
131 for (const auto i : indices_want)
132 {
134 {
135 indices_want_clean.push_back(static_cast<PetscInt>(i));
136 indices_want_loc.push_back(loc);
137 }
138 else
139 has_invalid = true;
140 loc++;
141 }
142 if (!has_invalid)
143 indices_want_loc.clear();
144
145 this->do_reinit(indices_has_clean,
146 indices_has_loc,
147 indices_want_clean,
148 indices_want_loc,
149 communicator);
150 }
151
152 void
153 CommunicationPattern::do_reinit(const std::vector<PetscInt> &inidx,
154 const std::vector<PetscInt> &inloc,
155 const std::vector<PetscInt> &outidx,
156 const std::vector<PetscInt> &outloc,
157 const MPI_Comm communicator)
158 {
159 clear();
160
161 // inidx is assumed to be unstructured and non-overlapping.
162 // However, it may have holes in it and not be a full cover.
163 //
164 // We create two PETSc SFs and compose them to get
165 // the final communication pattern
166 //
167 // sf1 : local distributed to tmp
168 // sf2 : tmp to local with ghosts
169 // sf(x) = sf2(sf1(x))
170 PetscSF sf1, sf2;
171
172 // First create an SF where leaves are inidx (at location inloc)
173 // and roots are unique indices in contiguous way
174 // Code adapted from MatZeroRowsMapLocal_Private in PETSc
175 PetscInt n = static_cast<PetscInt>(inidx.size());
176 PetscInt lN = n > 0 ? *std::max_element(inidx.begin(), inidx.end()) : -1;
177 PetscInt N, nl;
178
179 Utilities::MPI::internal::all_reduce<PetscInt>(
180 MPI_MAX,
182 communicator,
183 ArrayView<PetscInt>(&N, 1));
184
185 PetscSFNode *remotes;
186 AssertPETSc(PetscMalloc1(n, &remotes));
187
188 PetscLayout layout;
189 AssertPETSc(PetscLayoutCreate(communicator, &layout));
190 AssertPETSc(PetscLayoutSetSize(layout, N + 1));
191 AssertPETSc(PetscLayoutSetUp(layout));
192 AssertPETSc(PetscLayoutGetLocalSize(layout, &nl));
193
194 const PetscInt *ranges;
195 AssertPETSc(PetscLayoutGetRanges(layout, &ranges));
196
197 PetscInt cnt = 0;
198 PetscMPIInt owner = 0;
199 for (const auto idx : inidx)
200 {
201 // short-circuit the search if the last owner owns this index too
202 if (idx < ranges[owner] || ranges[owner + 1] <= idx)
203 {
204 AssertPETSc(PetscLayoutFindOwner(layout, idx, &owner));
205 }
206 remotes[cnt].rank = owner;
207 remotes[cnt].index = idx - ranges[owner];
208 cnt++;
209 }
210
211 AssertPETSc(PetscSFCreate(communicator, &sf2));
212 AssertPETSc(PetscSFSetGraph(sf2,
213 nl,
214 n,
215 const_cast<PetscInt *>(
216 inloc.size() > 0 ? inloc.data() : nullptr),
217 PETSC_COPY_VALUES,
218 remotes,
219 PETSC_OWN_POINTER));
220 AssertPETSc(PetscSFSetUp(sf2));
221 // We need to invert root and leaf space to create the first SF
222 AssertPETSc(PetscSFCreateInverseSF(sf2, &sf1));
223 AssertPETSc(PetscSFDestroy(&sf2));
224
225 // Now create the SF from the contiguous space to the local output space
226 n = static_cast<PetscInt>(outidx.size());
227 AssertPETSc(PetscSFCreate(communicator, &sf2));
228 AssertPETSc(PetscSFSetGraphLayout(
229 sf2,
230 layout,
231 n,
232 const_cast<PetscInt *>(outloc.size() > 0 ? outloc.data() : nullptr),
233 PETSC_COPY_VALUES,
234 const_cast<PetscInt *>(n > 0 ? outidx.data() : nullptr)));
235 AssertPETSc(PetscSFSetUp(sf2));
236
237 // The final SF is the composition of the two
238 AssertPETSc(PetscSFCompose(sf1, sf2, &sf));
239
240 // Cleanup
241 AssertPETSc(PetscLayoutDestroy(&layout));
242 AssertPETSc(PetscSFDestroy(&sf1));
243 AssertPETSc(PetscSFDestroy(&sf2));
244 }
245
246 void
248 {
249 AssertPETSc(PetscSFDestroy(&sf));
250 }
251
254 {
255 return PetscObjectComm(reinterpret_cast<PetscObject>(sf));
256 }
257
258 template <typename Number>
259 void
261 const ArrayView<const Number> &src,
262 const ArrayView<Number> & dst) const
263 {
264 auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
265
266# if DEAL_II_PETSC_VERSION_LT(3, 15, 0)
267 AssertPETSc(PetscSFBcastBegin(sf, datatype, src.data(), dst.data()));
268# else
270 PetscSFBcastBegin(sf, datatype, src.data(), dst.data(), MPI_REPLACE));
271# endif
272 }
273
274 template <typename Number>
275 void
277 const ArrayView<const Number> &src,
278 const ArrayView<Number> & dst) const
279 {
280 auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
281
282# if DEAL_II_PETSC_VERSION_LT(3, 15, 0)
283 AssertPETSc(PetscSFBcastEnd(sf, datatype, src.data(), dst.data()));
284# else
286 PetscSFBcastEnd(sf, datatype, src.data(), dst.data(), MPI_REPLACE));
287# endif
288 }
289
290 template <typename Number>
291 void
293 const ArrayView<const Number> &src,
294 const ArrayView<Number> & dst) const
295 {
298 }
299
300 template <typename Number>
301 void
304 const ArrayView<const Number> &src,
305 const ArrayView<Number> & dst) const
306 {
307 MPI_Op mpiop = (op == VectorOperation::insert) ? MPI_REPLACE : MPI_SUM;
308 auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
309
311 PetscSFReduceBegin(sf, datatype, src.data(), dst.data(), mpiop));
312 }
313
314 template <typename Number>
315 void
318 const ArrayView<const Number> &src,
319 const ArrayView<Number> & dst) const
320 {
321 MPI_Op mpiop = (op == VectorOperation::insert) ? MPI_REPLACE : MPI_SUM;
322 auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
323
324 AssertPETSc(PetscSFReduceEnd(sf, datatype, src.data(), dst.data(), mpiop));
325 }
326
327 template <typename Number>
328 void
331 const ArrayView<const Number> &src,
332 const ArrayView<Number> & dst) const
333 {
336 }
337
338 // Partitioner
339
341 : ghost()
342 , larger_ghost()
343 , ghost_indices_data()
344 , n_ghost_indices_data(numbers::invalid_dof_index)
345 , n_ghost_indices_larger(numbers::invalid_dof_index)
346 {}
347
348 void
349 Partitioner::reinit(const IndexSet &locally_owned_indices,
350 const IndexSet &ghost_indices,
351 const MPI_Comm communicator)
352 {
354 ghost_indices_data.subtract_set(locally_owned_indices);
356
357 ghost.reinit(locally_owned_indices, ghost_indices_data, communicator);
359
362 }
363
364 void
365 Partitioner::reinit(const IndexSet &locally_owned_indices,
366 const IndexSet &ghost_indices,
367 const IndexSet &larger_ghost_indices,
368 const MPI_Comm communicator)
369 {
370 std::vector<types::global_dof_index> local_indices;
371 locally_owned_indices.fill_index_vector(local_indices);
372
374 ghost_indices_data.subtract_set(locally_owned_indices);
376
377 std::vector<types::global_dof_index> expanded_ghost_indices(
378 larger_ghost_indices.n_elements(), numbers::invalid_dof_index);
379 for (auto index : ghost_indices_data)
380 {
381 Assert(larger_ghost_indices.is_element(index),
382 ExcMessage("The given larger ghost index set must contain "
383 "all indices in the actual index set."));
384 auto tmp_index = larger_ghost_indices.index_within_set(index);
385 expanded_ghost_indices[tmp_index] = index;
386 }
387
388 ghost.reinit(locally_owned_indices, ghost_indices_data, communicator);
389 larger_ghost.reinit(local_indices, expanded_ghost_indices, communicator);
391 n_ghost_indices_larger = larger_ghost_indices.n_elements();
392 }
393
396 {
398 }
399
400 template <typename Number>
401 void
403 const ArrayView<Number> &dst) const
404 {
405 if (dst.size() == n_ghost_indices_larger)
406 {
408 }
409 else
410 {
412 }
413 }
414
415 template <typename Number>
416 void
418 const ArrayView<const Number> &src,
419 const ArrayView<Number> & dst) const
420 {
421 if (dst.size() == n_ghost_indices_larger)
422 {
424 }
425 else
426 {
428 }
429 }
430
431 template <typename Number>
432 void
434 const ArrayView<Number> & dst) const
435 {
438 }
439
440 template <typename Number>
441 void
444 const ArrayView<const Number> &src,
445 const ArrayView<Number> & dst) const
446 {
447 if (src.size() == n_ghost_indices_larger)
448 {
450 }
451 else
452 {
454 }
455 }
456
457 template <typename Number>
458 void
461 const ArrayView<const Number> &src,
462 const ArrayView<Number> & dst) const
463 {
464 if (src.size() == n_ghost_indices_larger)
465 {
467 }
468 else
469 {
471 }
472 }
473
474 template <typename Number>
475 void
477 const ArrayView<const Number> &src,
478 const ArrayView<Number> &dst) const
479 {
482 }
483
484} // namespace PETScWrappers
485
486// Explicit instantiations
487# include "petsc_communication_pattern.inst"
488
490
491#endif // DEAL_II_WITH_PETSC
value_type * data() const noexcept
Definition array_view.h:553
std::size_t size() const
Definition array_view.h:576
IS make_petsc_is(const MPI_Comm communicator=MPI_COMM_WORLD) const
Definition index_set.cc:858
void fill_index_vector(std::vector< size_type > &indices) const
Definition index_set.cc:709
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1883
size_type n_elements() const
Definition index_set.h:1816
bool is_element(const size_type index) const
Definition index_set.h:1776
void subtract_set(const IndexSet &other)
Definition index_set.cc:268
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1697
void compress() const
Definition index_set.h:1669
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1727
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
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
virtual void reinit(const IndexSet &locally_owned_indices, const IndexSet &ghost_indices, const MPI_Comm communicator) override
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
const types::global_dof_index invalid_dof_index
Definition types.h:233
#define AssertPETSc(code)