Reference documentation for deal.II version GIT e8f57130d0 2023-09-28 07:00:03+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\}}\)
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 
22 # include <deal.II/lac/exceptions.h>
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 
34 namespace PETScWrappers
35 {
37  : sf(nullptr)
38  {}
39 
40 
41 
43  {
44  clear();
45  }
46 
47 
48 
49  void
51  const IndexSet &ghost_indices,
52  const MPI_Comm communicator)
53  {
54  clear();
55 
56  PetscLayout layout;
57  AssertPETSc(PetscLayoutCreate(communicator, &layout));
58  AssertPETSc(PetscLayoutSetLocalSize(layout, local_size));
59  AssertPETSc(PetscLayoutSetUp(layout));
60 
61  PetscInt start, end;
62  AssertPETSc(PetscLayoutGetRange(layout, &start, &end));
63 
64  IndexSet want;
65  want.add_range(start, end);
66  want.add_indices(ghost_indices);
67  want.compress();
68 
69  const PetscInt *idxs;
70  PetscInt n;
71  IS is = want.make_petsc_is(communicator);
72  AssertPETSc(ISGetLocalSize(is, &n));
73  AssertPETSc(ISGetIndices(is, &idxs));
74 
75  AssertPETSc(PetscSFCreate(communicator, &sf));
77  PetscSFSetGraphLayout(sf, layout, n, nullptr, PETSC_OWN_POINTER, idxs));
78  AssertPETSc(PetscSFSetUp(sf));
79 
80  AssertPETSc(ISRestoreIndices(is, &idxs));
81  AssertPETSc(ISDestroy(&is));
82  AssertPETSc(PetscLayoutDestroy(&layout));
83  }
84 
85 
86 
87  void
88  CommunicationPattern::reinit(const IndexSet &locally_owned_indices,
89  const IndexSet &ghost_indices,
90  const MPI_Comm communicator)
91  {
92  std::vector<types::global_dof_index> in_deal;
93  locally_owned_indices.fill_index_vector(in_deal);
94  std::vector<PetscInt> in_petsc(in_deal.begin(), in_deal.end());
95 
96  std::vector<types::global_dof_index> out_deal;
97  ghost_indices.fill_index_vector(out_deal);
98  std::vector<PetscInt> out_petsc(out_deal.begin(), out_deal.end());
99 
100  std::vector<PetscInt> dummy;
101 
102  this->do_reinit(in_petsc, dummy, out_petsc, dummy, communicator);
103  }
104 
105 
106 
107  void
109  const std::vector<types::global_dof_index> &indices_has,
110  const std::vector<types::global_dof_index> &indices_want,
111  const MPI_Comm communicator)
112  {
113  // Clean vectors from numbers::invalid_dof_index (indicating padding)
114  std::vector<PetscInt> indices_has_clean, indices_has_loc;
115  std::vector<PetscInt> indices_want_clean, indices_want_loc;
116  indices_want_clean.reserve(indices_want.size());
117  indices_want_loc.reserve(indices_want.size());
118  indices_has_clean.reserve(indices_has.size());
119  indices_has_loc.reserve(indices_has.size());
120 
121  PetscInt loc = 0;
122  bool has_invalid = false;
123  for (const auto i : indices_has)
124  {
126  {
127  indices_has_clean.push_back(static_cast<PetscInt>(i));
128  indices_has_loc.push_back(loc);
129  }
130  else
131  has_invalid = true;
132  loc++;
133  }
134  if (!has_invalid)
135  indices_has_loc.clear();
136 
137  loc = 0;
138  has_invalid = false;
139  for (const auto i : indices_want)
140  {
142  {
143  indices_want_clean.push_back(static_cast<PetscInt>(i));
144  indices_want_loc.push_back(loc);
145  }
146  else
147  has_invalid = true;
148  loc++;
149  }
150  if (!has_invalid)
151  indices_want_loc.clear();
152 
153  this->do_reinit(indices_has_clean,
154  indices_has_loc,
155  indices_want_clean,
156  indices_want_loc,
157  communicator);
158  }
159 
160 
161 
162  void
163  CommunicationPattern::do_reinit(const std::vector<PetscInt> &inidx,
164  const std::vector<PetscInt> &inloc,
165  const std::vector<PetscInt> &outidx,
166  const std::vector<PetscInt> &outloc,
167  const MPI_Comm communicator)
168  {
169  clear();
170 
171  // inidx is assumed to be unstructured and non-overlapping.
172  // However, it may have holes in it and not be a full cover.
173  //
174  // We create two PETSc SFs and compose them to get
175  // the final communication pattern
176  //
177  // sf1 : local distributed to tmp
178  // sf2 : tmp to local with ghosts
179  // sf(x) = sf2(sf1(x))
180  PetscSF sf1, sf2;
181 
182  // First create an SF where leaves are inidx (at location inloc)
183  // and roots are unique indices in contiguous way
184  // Code adapted from MatZeroRowsMapLocal_Private in PETSc
185  PetscInt n = static_cast<PetscInt>(inidx.size());
186  PetscInt lN = n > 0 ? *std::max_element(inidx.begin(), inidx.end()) : -1;
187  PetscInt N, nl;
188 
189  Utilities::MPI::internal::all_reduce<PetscInt>(
190  MPI_MAX,
192  communicator,
193  ArrayView<PetscInt>(&N, 1));
194 
195  PetscSFNode *remotes;
196  AssertPETSc(PetscMalloc1(n, &remotes));
197 
198  PetscLayout layout;
199  AssertPETSc(PetscLayoutCreate(communicator, &layout));
200  AssertPETSc(PetscLayoutSetSize(layout, N + 1));
201  AssertPETSc(PetscLayoutSetUp(layout));
202  AssertPETSc(PetscLayoutGetLocalSize(layout, &nl));
203 
204  const PetscInt *ranges;
205  AssertPETSc(PetscLayoutGetRanges(layout, &ranges));
206 
207  PetscInt cnt = 0;
208  PetscMPIInt owner = 0;
209  for (const auto idx : inidx)
210  {
211  // short-circuit the search if the last owner owns this index too
212  if (idx < ranges[owner] || ranges[owner + 1] <= idx)
213  {
214  AssertPETSc(PetscLayoutFindOwner(layout, idx, &owner));
215  }
216  remotes[cnt].rank = owner;
217  remotes[cnt].index = idx - ranges[owner];
218  cnt++;
219  }
220 
221  AssertPETSc(PetscSFCreate(communicator, &sf2));
222  AssertPETSc(PetscSFSetGraph(sf2,
223  nl,
224  n,
225  const_cast<PetscInt *>(
226  inloc.size() > 0 ? inloc.data() : nullptr),
227  PETSC_COPY_VALUES,
228  remotes,
229  PETSC_OWN_POINTER));
230  AssertPETSc(PetscSFSetUp(sf2));
231  // We need to invert root and leaf space to create the first SF
232  AssertPETSc(PetscSFCreateInverseSF(sf2, &sf1));
233  AssertPETSc(PetscSFDestroy(&sf2));
234 
235  // Now create the SF from the contiguous space to the local output space
236  n = static_cast<PetscInt>(outidx.size());
237  AssertPETSc(PetscSFCreate(communicator, &sf2));
238  AssertPETSc(PetscSFSetGraphLayout(
239  sf2,
240  layout,
241  n,
242  const_cast<PetscInt *>(outloc.size() > 0 ? outloc.data() : nullptr),
243  PETSC_COPY_VALUES,
244  const_cast<PetscInt *>(n > 0 ? outidx.data() : nullptr)));
245  AssertPETSc(PetscSFSetUp(sf2));
246 
247  // The final SF is the composition of the two
248  AssertPETSc(PetscSFCompose(sf1, sf2, &sf));
249 
250  // Cleanup
251  AssertPETSc(PetscLayoutDestroy(&layout));
252  AssertPETSc(PetscSFDestroy(&sf1));
253  AssertPETSc(PetscSFDestroy(&sf2));
254  }
255 
256 
257 
258  void
260  {
261  AssertPETSc(PetscSFDestroy(&sf));
262  }
263 
264 
265 
266  MPI_Comm
268  {
269  return PetscObjectComm(reinterpret_cast<PetscObject>(sf));
270  }
271 
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(PetscSFBcastBegin(sf, datatype, src.data(), dst.data()));
284 # else
285  AssertPETSc(
286  PetscSFBcastBegin(sf, datatype, src.data(), dst.data(), MPI_REPLACE));
287 # endif
288  }
289 
290 
291 
292  template <typename Number>
293  void
295  const ArrayView<const Number> &src,
296  const ArrayView<Number> &dst) const
297  {
298  auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
299 
300 # if DEAL_II_PETSC_VERSION_LT(3, 15, 0)
301  AssertPETSc(PetscSFBcastEnd(sf, datatype, src.data(), dst.data()));
302 # else
303  AssertPETSc(
304  PetscSFBcastEnd(sf, datatype, src.data(), dst.data(), MPI_REPLACE));
305 # endif
306  }
307 
308 
309 
310  template <typename Number>
311  void
313  const ArrayView<const Number> &src,
314  const ArrayView<Number> &dst) const
315  {
318  }
319 
320 
321 
322  template <typename Number>
323  void
325  const VectorOperation::values op,
326  const ArrayView<const Number> &src,
327  const ArrayView<Number> &dst) const
328  {
329  MPI_Op mpiop = (op == VectorOperation::insert) ? MPI_REPLACE : MPI_SUM;
330  auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
331 
332  AssertPETSc(
333  PetscSFReduceBegin(sf, datatype, src.data(), dst.data(), mpiop));
334  }
335 
336 
337 
338  template <typename Number>
339  void
341  const VectorOperation::values op,
342  const ArrayView<const Number> &src,
343  const ArrayView<Number> &dst) const
344  {
345  MPI_Op mpiop = (op == VectorOperation::insert) ? MPI_REPLACE : MPI_SUM;
346  auto datatype = Utilities::MPI::mpi_type_id_for_type<Number>;
347 
348  AssertPETSc(PetscSFReduceEnd(sf, datatype, src.data(), dst.data(), mpiop));
349  }
350 
351 
352 
353  template <typename Number>
354  void
356  const VectorOperation::values op,
357  const ArrayView<const Number> &src,
358  const ArrayView<Number> &dst) const
359  {
360  import_from_ghosted_array_start(op, src, dst);
361  import_from_ghosted_array_finish(op, src, dst);
362  }
363 
364 
365 
366  // Partitioner
367 
369  : ghost()
370  , larger_ghost()
371  , ghost_indices_data()
372  , n_ghost_indices_data(numbers::invalid_dof_index)
373  , n_ghost_indices_larger(numbers::invalid_dof_index)
374  {}
375 
376 
377 
378  void
379  Partitioner::reinit(const IndexSet &locally_owned_indices,
380  const IndexSet &ghost_indices,
381  const MPI_Comm communicator)
382  {
384  ghost_indices_data.subtract_set(locally_owned_indices);
386 
387  ghost.reinit(locally_owned_indices, ghost_indices_data, communicator);
389 
392  }
393 
394  void
395  Partitioner::reinit(const IndexSet &locally_owned_indices,
396  const IndexSet &ghost_indices,
397  const IndexSet &larger_ghost_indices,
398  const MPI_Comm communicator)
399  {
400  std::vector<types::global_dof_index> local_indices;
401  locally_owned_indices.fill_index_vector(local_indices);
402 
404  ghost_indices_data.subtract_set(locally_owned_indices);
406 
407  std::vector<types::global_dof_index> expanded_ghost_indices(
408  larger_ghost_indices.n_elements(), numbers::invalid_dof_index);
409  for (auto index : ghost_indices_data)
410  {
411  Assert(larger_ghost_indices.is_element(index),
412  ExcMessage("The given larger ghost index set must contain "
413  "all indices in the actual index set."));
414  auto tmp_index = larger_ghost_indices.index_within_set(index);
415  expanded_ghost_indices[tmp_index] = index;
416  }
417 
418  ghost.reinit(locally_owned_indices, ghost_indices_data, communicator);
419  larger_ghost.reinit(local_indices, expanded_ghost_indices, communicator);
421  n_ghost_indices_larger = larger_ghost_indices.n_elements();
422  }
423 
424  MPI_Comm
426  {
427  return ghost.get_mpi_communicator();
428  }
429 
430  template <typename Number>
431  void
433  const ArrayView<Number> &dst) const
434  {
435  if (dst.size() == n_ghost_indices_larger)
436  {
438  }
439  else
440  {
442  }
443  }
444 
445  template <typename Number>
446  void
448  const ArrayView<const Number> &src,
449  const ArrayView<Number> &dst) const
450  {
451  if (dst.size() == n_ghost_indices_larger)
452  {
454  }
455  else
456  {
458  }
459  }
460 
461  template <typename Number>
462  void
464  const ArrayView<Number> &dst) const
465  {
468  }
469 
470  template <typename Number>
471  void
473  const VectorOperation::values op,
474  const ArrayView<const Number> &src,
475  const ArrayView<Number> &dst) const
476  {
477  if (src.size() == n_ghost_indices_larger)
478  {
480  }
481  else
482  {
484  }
485  }
486 
487  template <typename Number>
488  void
490  const VectorOperation::values op,
491  const ArrayView<const Number> &src,
492  const ArrayView<Number> &dst) const
493  {
494  if (src.size() == n_ghost_indices_larger)
495  {
497  }
498  else
499  {
501  }
502  }
503 
504  template <typename Number>
505  void
507  const ArrayView<const Number> &src,
508  const ArrayView<Number> &dst) const
509  {
510  import_from_ghosted_array_start(op, src, dst);
511  import_from_ghosted_array_finish(op, src, dst);
512  }
513 
514 } // namespace PETScWrappers
515 
516 // Explicit instantiations
517 # include "petsc_communication_pattern.inst"
518 
520 
521 #endif // DEAL_II_WITH_PETSC
value_type * data() const noexcept
Definition: array_view.h:550
std::size_t size() const
Definition: array_view.h:573
IS make_petsc_is(const MPI_Comm communicator=MPI_COMM_WORLD) const
Definition: index_set.cc:909
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:760
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
size_type n_elements() const
Definition: index_set.h:1854
bool is_element(const size_type index) const
Definition: index_set.h:1814
void subtract_set(const IndexSet &other)
Definition: index_set.cc:288
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1723
void compress() const
Definition: index_set.h:1704
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1751
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
const IndexSet & ghost_indices() const
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcMessage(std::string arg1)
static const char N
VectorType::value_type * end(VectorType &V)
const types::global_dof_index invalid_dof_index
Definition: types.h:233
unsigned int global_dof_index
Definition: types.h:82
#define AssertPETSc(code)