Reference documentation for deal.II version GIT 7b2de2f2f9 2023-09-24 11:00:02+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_parallel_vector.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 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 
16 #include <deal.II/base/mpi.h>
17 
19 
20 #ifdef DEAL_II_WITH_PETSC
21 
22 # include <algorithm>
23 # include <cmath>
24 
26 
27 namespace PETScWrappers
28 {
29  namespace MPI
30  {
31  // Check that the class we declare here satisfies the
32  // vector-space-vector concept. If we catch it here,
33  // any mistake in the vector class declaration would
34  // show up in uses of this class later on as well.
35 # ifdef DEAL_II_HAVE_CXX20
36  static_assert(concepts::is_vector_space_vector<Vector>);
37 # endif
38 
40  {
41  // virtual functions called in constructors and destructors never use the
42  // override in a derived class
43  // for clarity be explicit on which function is called
44  Vector::create_vector(MPI_COMM_SELF, 0, 0);
45  }
46 
47 
48 
49  Vector::Vector(const MPI_Comm communicator,
50  const size_type n,
51  const size_type locally_owned_size)
52  {
53  Vector::create_vector(communicator, n, locally_owned_size);
54  }
55 
56 
57 
58  Vector::Vector(const IndexSet &local,
59  const IndexSet &ghost,
60  const MPI_Comm communicator)
61  {
62  Assert(local.is_ascending_and_one_to_one(communicator),
64 
65  IndexSet ghost_set = ghost;
66  ghost_set.subtract_set(local);
67 
68  Vector::create_vector(communicator,
69  local.size(),
70  local.n_elements(),
71  ghost_set);
72  }
73 
74 
75 
77  : VectorBase()
78  {
79  if (v.has_ghost_elements())
81  v.size(),
83  v.ghost_indices);
84  else
86  v.size(),
87  v.locally_owned_size());
88 
89  this->operator=(v);
90  }
91 
92 
93 
94  Vector::Vector(const IndexSet &local, const MPI_Comm communicator)
95  {
96  Assert(local.is_ascending_and_one_to_one(communicator),
98  Vector::create_vector(communicator, local.size(), local.n_elements());
99  }
100 
101 
102 
103  Vector &
105  {
106  // make sure left- and right-hand side of the assignment are
107  // compress()'ed:
109  internal::VectorReference::ExcWrongMode(VectorOperation::unknown,
110  v.last_action));
112  internal::VectorReference::ExcWrongMode(VectorOperation::unknown,
113  last_action));
114 
115  // if the vectors have different sizes,
116  // then first resize the present one
117  if (size() != v.size())
118  {
119  if (v.has_ghost_elements())
121  v.ghost_indices,
123  else
125  v.size(),
126  v.locally_owned_size(),
127  true);
128  }
129 
130  PetscErrorCode ierr = VecCopy(v.vector, vector);
131  AssertThrow(ierr == 0, ExcPETScError(ierr));
132 
133  if (has_ghost_elements())
134  {
135  ierr = VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
136  AssertThrow(ierr == 0, ExcPETScError(ierr));
137  ierr = VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
138  AssertThrow(ierr == 0, ExcPETScError(ierr));
139  }
140  return *this;
141  }
142 
143 
144 
145  void
147  {
149 
150  create_vector(MPI_COMM_SELF, 0, 0);
151  }
152 
153 
154 
155  void
156  Vector::reinit(const MPI_Comm communicator,
157  const size_type n,
158  const size_type local_sz,
159  const bool omit_zeroing_entries)
160  {
161  // only do something if the sizes
162  // mismatch (may not be true for every proc)
163 
164  int k_global, k = ((size() != n) || (locally_owned_size() != local_sz));
165  {
166  const int ierr =
167  MPI_Allreduce(&k, &k_global, 1, MPI_INT, MPI_LOR, communicator);
168  AssertThrowMPI(ierr);
169  }
170 
171  if (k_global || has_ghost_elements())
172  {
173  // FIXME: I'd like to use this here,
174  // but somehow it leads to odd errors
175  // somewhere down the line in some of
176  // the tests:
177  // const PetscErrorCode ierr = VecSetSizes (vector, n, n);
178  // AssertThrow (ierr == 0, ExcPETScError(ierr));
179 
180  // so let's go the slow way:
181 
182  const PetscErrorCode ierr = VecDestroy(&vector);
183  AssertThrow(ierr == 0, ExcPETScError(ierr));
184 
185  create_vector(communicator, n, local_sz);
186  }
187 
188  // finally clear the new vector if so
189  // desired
190  if (omit_zeroing_entries == false)
191  *this = 0;
192  }
193 
194 
195 
196  void
197  Vector::reinit(const Vector &v, const bool omit_zeroing_entries)
198  {
199  if (v.has_ghost_elements())
200  {
202  v.ghost_indices,
204  if (!omit_zeroing_entries)
205  {
206  const PetscErrorCode ierr = VecSet(vector, 0.0);
207  AssertThrow(ierr == 0, ExcPETScError(ierr));
208  }
209  }
210  else
212  v.size(),
213  v.locally_owned_size(),
214  omit_zeroing_entries);
215  }
216 
217 
218 
219  void
220  Vector::reinit(const IndexSet &local,
221  const IndexSet &ghost,
222  const MPI_Comm comm)
223  {
224  const PetscErrorCode ierr = VecDestroy(&vector);
225  AssertThrow(ierr == 0, ExcPETScError(ierr));
226 
228 
229  IndexSet ghost_set = ghost;
230  ghost_set.subtract_set(local);
231 
232  create_vector(comm, local.size(), local.n_elements(), ghost_set);
233  }
234 
235  void
236  Vector::reinit(const IndexSet &local, const MPI_Comm comm)
237  {
238  const PetscErrorCode ierr = VecDestroy(&vector);
239  AssertThrow(ierr == 0, ExcPETScError(ierr));
240 
242  Assert(local.size() > 0, ExcMessage("can not create vector of size 0."));
243  create_vector(comm, local.size(), local.n_elements());
244  }
245 
246  void
248  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
249  const bool make_ghosted)
250  {
251  if (make_ghosted)
252  {
253  Assert(partitioner->ghost_indices_initialized(),
254  ExcMessage("You asked to create a ghosted vector, but the "
255  "partitioner does not provide ghost indices."));
256 
257  this->reinit(partitioner->locally_owned_range(),
258  partitioner->ghost_indices(),
259  partitioner->get_mpi_communicator());
260  }
261  else
262  {
263  this->reinit(partitioner->locally_owned_range(),
264  partitioner->get_mpi_communicator());
265  }
266  }
267 
268 
269  void
270  Vector::create_vector(const MPI_Comm communicator,
271  const size_type n,
272  const size_type locally_owned_size)
273  {
274  (void)n;
276  ghosted = false;
277 
278  const PetscErrorCode ierr = VecCreateMPI(communicator,
280  PETSC_DETERMINE,
281  &vector);
282  AssertThrow(ierr == 0, ExcPETScError(ierr));
283 
284  Assert(size() == n, ExcDimensionMismatch(size(), n));
285  }
286 
287 
288 
289  void
290  Vector::create_vector(const MPI_Comm communicator,
291  const size_type n,
292  const size_type locally_owned_size,
293  const IndexSet &ghostnodes)
294  {
295  (void)n;
297  ghosted = true;
298  ghost_indices = ghostnodes;
299 
300  const std::vector<size_type> ghostindices = ghostnodes.get_index_vector();
301 
302  const PetscInt *ptr =
303  (ghostindices.size() > 0 ?
304  reinterpret_cast<const PetscInt *>(ghostindices.data()) :
305  nullptr);
306 
307  PetscErrorCode ierr = VecCreateGhost(communicator,
309  PETSC_DETERMINE,
310  ghostindices.size(),
311  ptr,
312  &vector);
313  AssertThrow(ierr == 0, ExcPETScError(ierr));
314 
315  Assert(size() == n, ExcDimensionMismatch(size(), n));
316 
317 # ifdef DEBUG
318  {
319  // test ghost allocation in debug mode
320  PetscInt begin, end;
321 
322  ierr = VecGetOwnershipRange(vector, &begin, &end);
323  AssertThrow(ierr == 0, ExcPETScError(ierr));
324 
326  static_cast<size_type>(end - begin));
327 
328  Vec l;
329  ierr = VecGhostGetLocalForm(vector, &l);
330  AssertThrow(ierr == 0, ExcPETScError(ierr));
331 
332  PetscInt lsize;
333  ierr = VecGetSize(l, &lsize);
334  AssertThrow(ierr == 0, ExcPETScError(ierr));
335 
336  ierr = VecGhostRestoreLocalForm(vector, &l);
337  AssertThrow(ierr == 0, ExcPETScError(ierr));
338 
339  AssertDimension(lsize,
340  end - begin +
341  static_cast<PetscInt>(ghost_indices.n_elements()));
342  }
343 # endif
344  }
345 
346 
347 
348  bool
350  {
351  unsigned int has_nonzero = VectorBase::all_zero() ? 0 : 1;
352  // in parallel, check that the vector
353  // is zero on _all_ processors.
354  unsigned int num_nonzero =
355  Utilities::MPI::sum(has_nonzero, this->get_mpi_communicator());
356  return num_nonzero == 0;
357  }
358 
359 
360  void
361  Vector::print(std::ostream &out,
362  const unsigned int precision,
363  const bool scientific,
364  const bool across) const
365  {
366  AssertThrow(out.fail() == false, ExcIO());
367 
368  // get a representation of the vector and
369  // loop over all the elements
370  const PetscScalar *val;
371  PetscInt nlocal, istart, iend;
372 
373  PetscErrorCode ierr = VecGetArrayRead(vector, &val);
374  AssertThrow(ierr == 0, ExcPETScError(ierr));
375 
376  ierr = VecGetLocalSize(vector, &nlocal);
377  AssertThrow(ierr == 0, ExcPETScError(ierr));
378 
379  ierr = VecGetOwnershipRange(vector, &istart, &iend);
380  AssertThrow(ierr == 0, ExcPETScError(ierr));
381 
382  // save the state of out stream
383  std::ios::fmtflags old_flags = out.flags();
384  unsigned int old_precision = out.precision(precision);
385 
386  out.precision(precision);
387  if (scientific)
388  out.setf(std::ios::scientific, std::ios::floatfield);
389  else
390  out.setf(std::ios::fixed, std::ios::floatfield);
391 
392  // let each processor produce its output in turn. this requires
393  // synchronizing output between processors using a barrier --
394  // which is clearly slow, but nobody is going to print a whole
395  // matrix this way on a regular basis for production runs, so
396  // the slowness of the barrier doesn't matter
397  MPI_Comm communicator = this->get_mpi_communicator();
398  for (unsigned int i = 0;
399  i < Utilities::MPI::n_mpi_processes(communicator);
400  i++)
401  {
402  const int mpi_ierr = MPI_Barrier(communicator);
403  AssertThrowMPI(mpi_ierr);
404 
405  if (i == Utilities::MPI::this_mpi_process(communicator))
406  {
407  if (across)
408  {
409  out << "[Proc" << i << " " << istart << "-" << iend - 1 << "]"
410  << ' ';
411  for (PetscInt i = 0; i < nlocal; ++i)
412  out << val[i] << ' ';
413  }
414  else
415  {
416  out << "[Proc " << i << " " << istart << "-" << iend - 1
417  << "]" << std::endl;
418  for (PetscInt i = 0; i < nlocal; ++i)
419  out << val[i] << std::endl;
420  }
421  out << std::endl;
422  }
423  }
424  // reset output format
425  out.flags(old_flags);
426  out.precision(old_precision);
427 
428  // restore the representation of the
429  // vector
430  ierr = VecRestoreArrayRead(vector, &val);
431  AssertThrow(ierr == 0, ExcPETScError(ierr));
432 
433  AssertThrow(out.fail() == false, ExcIO());
434  }
435 
436  } // namespace MPI
437 
438 } // namespace PETScWrappers
439 
441 
442 #endif // DEAL_II_WITH_PETSC
bool is_ascending_and_one_to_one(const MPI_Comm communicator) const
Definition: index_set.cc:922
size_type size() const
Definition: index_set.h:1696
size_type n_elements() const
Definition: index_set.h:1854
void subtract_set(const IndexSet &other)
Definition: index_set.cc:288
std::vector< size_type > get_index_vector() const
Definition: index_set.cc:734
types::global_dof_index size_type
Definition: petsc_vector.h:164
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
virtual void clear() override
Vector & operator=(const Vector &v)
virtual void create_vector(const MPI_Comm comm, const size_type n, const size_type locally_owned_size)
VectorOperation::values last_action
IndexSet locally_owned_elements() const
MPI_Comm get_mpi_communicator() const
bool has_ghost_elements() const
size_type locally_owned_size() const
size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1916
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition: mpi.cc:149
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition: mpi.cc:164
const MPI_Comm comm