AMDiS  0.3
The Adaptive Multi-Dimensional Simulation Toolbox
VectorBackend.hpp
1 #pragma once
2 
3 #include <petscvec.h>
4 
5 #include <dune/common/ftraits.hh>
6 #include <dune/common/rangeutilities.hh>
7 #include <dune/common/shared_ptr.hh>
8 #include <dune/functions/functionspacebases/flatmultiindex.hh>
9 
10 #include <amdis/Output.hpp>
11 #include <amdis/common/FakeContainer.hpp>
12 #include <amdis/operations/Assigner.hpp>
13 #include <amdis/typetree/MultiIndex.hpp>
14 
15 namespace AMDiS
16 {
18  template <class DofMap>
20  {
21  template <class A>
22  using AssignMode = std::conditional_t<std::is_same_v<A,Assigner::plus_assign>,
23  std::integral_constant<::InsertMode, ADD_VALUES>,
24  std::integral_constant<::InsertMode, INSERT_VALUES>>;
25 
26  public:
28  using BaseVector = ::Vec;
29 
31  using value_type = PetscScalar;
32 
34  using size_type = PetscInt;
35 
36  public:
38  template <class Basis>
39  explicit PetscVector(Basis const& basis)
40  : dofMap_(basis.communicator().dofMap())
41  {}
42 
45  : dofMap_(other.dofMap_)
46  {
47  swap(*this, other);
48  }
49 
51  PetscVector(PetscVector const& other)
52  : dofMap_(other.dofMap_)
53  , initialized_(other.initialized_)
54  {
55  if (other.initialized_) {
56  VecDuplicate(other.vector_, &vector_);
57  VecCopy(other.vector_, vector_);
58  }
59  }
60 
63  {
64  destroy();
65  }
66 
69  {
70  swap(*this, other);
71  return *this;
72  }
73 
76  {
77  return *this = PetscVector(other);
78  }
79 
81  BaseVector const& vector() const
82  {
83  return vector_;
84  }
85 
88  {
89  return vector_;
90  }
91 
93  std::size_t localSize() const
94  {
95  return dofMap_.localSize();
96  }
97 
99  std::size_t globalSize() const
100  {
101  return dofMap_.globalSize();
102  }
103 
105  // [[possibly collective]]
106  template <class SizeInfo>
107  void init(SizeInfo const&, bool clear)
108  {
109  Dune::Timer t;
110  if (!initialized_) {
111  // create the vector
112  create();
113  info(3, " create new vector needed {} seconds", t.elapsed());
114  t.reset();
115  }
116  else {
117  PetscInt localSize = 0, globalSize = 0;
118  VecGetSize(vector_, &globalSize);
119  VecGetLocalSize(vector_, &localSize);
120 
121  // re-create the vector
122  if (std::size_t(localSize) != dofMap_.localSize() ||
123  std::size_t(globalSize) != dofMap_.globalSize()) {
124  destroy();
125  create();
126  info(3, " destroy old and create new vector needed {} seconds", t.elapsed());
127  t.reset();
128  }
129  }
130 
131  if (clear)
132  setZero();
133 
134  initialized_ = true;
135  }
136 
138  // [[collective]]
139  void finish()
140  {
141  Dune::Timer t;
142  VecAssemblyBegin(vector_);
143  VecAssemblyEnd(vector_);
144  info(3, " finish vector assembly needed {} seconds", t.elapsed());
145  }
146 
148  // [[collective]]
149  void synchronize()
150  {
151  Dune::Timer t;
152  VecGhostUpdateBegin(vector_, INSERT_VALUES, SCATTER_FORWARD);
153  VecGhostUpdateEnd(vector_, INSERT_VALUES, SCATTER_FORWARD);
154  info(3, " synchronize vector needed {} seconds", t.elapsed());
155  }
156 
158  // [[not collective]]
159  template <class MultiIndex>
160  PetscScalar at(MultiIndex const& idx) const
161  {
162  PetscScalar y = 0;
163  gather(std::array<MultiIndex,1>{idx}, &y);
164  return y;
165  }
166 
168  // [[not collective]]
169  template <class MultiIndex, class Assign>
170  void insert(MultiIndex const& dof, PetscScalar value, Assign)
171  {
172  assert(initialized_);
173  VecSetValue(vector_, dofMap_.global(flatMultiIndex(dof)), value, AssignMode<Assign>::value);
174  }
175 
177  // [[not collective]]
178  template <class IndexRange, class OutputIterator>
179  void gather(IndexRange const& localInd, OutputIterator buffer) const
180  {
181  forEach(localInd, [&buffer](auto&&, auto const& value) { *buffer = value; ++buffer; });
182  }
183 
185 
189  // [[not collective]]
190  template <class IndexRange, class LocalValues, class MaskRange, class Assign>
191  void scatter(IndexRange const& localInd, LocalValues const& localVal, MaskRange const& mask, Assign)
192  {
193  static_assert(std::is_same_v<typename LocalValues::value_type, PetscScalar>, "");
194  assert(initialized_);
195  assert(localInd.size() == localVal.size());
196 
197  // map local to global indices
198  std::vector<PetscInt> globalInd;
199  globalInd.reserve(localInd.size());
200 
201  auto ind_it = std::begin(localInd);
202  auto mask_it = std::begin(mask);
203  for (; ind_it != std::end(localInd); ++mask_it, ++ind_it)
204  globalInd.push_back(*mask_it ? dofMap_.global(flatMultiIndex(*ind_it)) : -1);
205 
206  if (!std::is_same_v<MaskRange, FakeContainer<bool,true>>)
207  VecSetOption(vector_, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
208 
209  VecSetValues(vector_, globalInd.size(), globalInd.data(), localVal.data(), AssignMode<Assign>::value);
210  VecSetOption(vector_, VEC_IGNORE_NEGATIVE_INDICES, PETSC_FALSE);
211  }
212 
214  // [[not collective]]
215  template <class IndexRange, class Func>
216  void forEach(IndexRange const& localInd, Func&& f) const
217  {
218  // Obtains the local ghosted representation of a parallel vector
219  Vec localGhostedVec = nullptr;
220  VecGhostGetLocalForm(vector_, &localGhostedVec);
221  assert(localGhostedVec != nullptr);
222 
223  // A pointer to a contiguous array that contains this processor's portion of the vector data.
224  PetscScalar *ptr;
225  VecGetArray(localGhostedVec, &ptr);
226 
227  for (auto const& dof : localInd)
228  {
229  size_type i = flatMultiIndex(dof);
230  if (dofMap_.owner(i)) {
231  // the index is a purely local entry
232  f(dof, ptr[dofMap_.dofToLocal(i)]);
233  }
234  else {
235  // ghost entries are stored at the end of the array
236  f(dof, ptr[dofMap_.localSize() + dofMap_.dofToGhost(i)]);
237  }
238  }
239 
240  VecRestoreArray(localGhostedVec, &ptr);
241  VecGhostRestoreLocalForm(vector_, &localGhostedVec);
242  }
243 
244  template <class Func>
245  void forEach(Func&& f)
246  {
247  forEach(mappedRangeView(Dune::range(localSize()),
248  [](auto i) { return Dune::Functions::FlatMultiIndex<decltype(i)>{i}; }), FWD(f));
249  }
250 
251  template <class Func>
252  void forEach(Func&& f) const
253  {
254  forEach(mappedRangeView(Dune::range(localSize()),
255  [](auto i) { return Dune::Functions::FlatMultiIndex<decltype(i)>{i}; }), FWD(f));
256  }
257 
258 
260  // [[not collective]]
261  void set(PetscScalar value)
262  {
263  Vec localGhostedVec = nullptr;
264  VecGhostGetLocalForm(vector_, &localGhostedVec);
265  assert(localGhostedVec != nullptr);
266 
267  VecSet(localGhostedVec, value);
268  VecGhostRestoreLocalForm(vector_, &localGhostedVec);
269  }
270 
272  // [[not collective]]
273  void setZero()
274  {
275  Vec localGhostedVec = nullptr;
276  VecGhostGetLocalForm(vector_, &localGhostedVec);
277  assert(localGhostedVec != nullptr);
278 
279  VecZeroEntries(localGhostedVec);
280  VecGhostRestoreLocalForm(vector_, &localGhostedVec);
281  }
282 
283  friend void swap(PetscVector& lhs, PetscVector& rhs)
284  {
285  assert(&lhs.dofMap_ == &rhs.dofMap_);
286  std::swap(lhs.vector_, rhs.vector_);
287  std::swap(lhs.initialized_, rhs.initialized_);
288  }
289 
290  private:
291  // An MPI Communicator or PETSC_COMM_SELF
292  MPI_Comm comm() const
293  {
294 #if HAVE_MPI
295  return Environment::comm();
296 #else
297  return PETSC_COMM_SELF;
298 #endif
299  }
300 
301  // Create a new Vector
302  void create()
303  {
304  VecCreateGhost(comm(),
305  dofMap_.localSize(), dofMap_.globalSize(), dofMap_.ghostSize(), dofMap_.ghostIndices().data(),
306  &vector_);
307  }
308 
309  // Destroy an old vector if created before
310  void destroy()
311  {
312  if (initialized_)
313  VecDestroy(&vector_);
314 
315  initialized_ = false;
316  }
317 
318  private:
319  // The local-to-global map
320  DofMap const& dofMap_;
321 
323  mutable Vec vector_ = nullptr;
324  bool initialized_ = false;
325  };
326 
327 } // end namespace AMDiS
constexpr bool MultiIndex
A multi-index type.
Definition: Concepts.hpp:150
BaseVector const & vector() const
Return the data-vector vector_.
Definition: VectorBackend.hpp:81
PetscVector(PetscVector &&other)
Move constructor.
Definition: VectorBackend.hpp:44
BaseVector & vector()
Return the data-vector vector_.
Definition: VectorBackend.hpp:87
PetscVector(PetscVector const &other)
Copy constructor.
Definition: VectorBackend.hpp:51
~PetscVector()
Destructor.
Definition: VectorBackend.hpp:62
std::size_t globalSize() const
Return the number of entries in the global vector.
Definition: VectorBackend.hpp:99
void forEach(IndexRange const &localInd, Func &&f) const
Apply the functor f to each vector entry to local index in localInd range.
Definition: VectorBackend.hpp:216
PetscScalar value_type
The type of the elements of the DOFVector.
Definition: VectorBackend.hpp:31
void insert(MultiIndex const &dof, PetscScalar value, Assign)
Access the entry i of the vector with write-access.
Definition: VectorBackend.hpp:170
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
PetscVector & operator=(PetscVector &&other)
Move assignment operator.
Definition: VectorBackend.hpp:68
auto mappedRangeView(R &&range, F const &f)
Create a MappedRangeView.
Definition: MappedRangeView.hpp:445
PetscVector & operator=(PetscVector const &other)
Copy assignment operator.
Definition: VectorBackend.hpp:75
void init(SizeInfo const &, bool clear)
Resize the vector_ to the size s.
Definition: VectorBackend.hpp:107
void setZero()
Zero all entries, including the ghost entries.
Definition: VectorBackend.hpp:273
void scatter(IndexRange const &localInd, LocalValues const &localVal, MaskRange const &mask, Assign)
Store all values given by localVal with indices localInd in the vector.
Definition: VectorBackend.hpp:191
std::size_t localSize() const
Return the number of entries in the local part of the vector.
Definition: VectorBackend.hpp:93
void gather(IndexRange const &localInd, OutputIterator buffer) const
Collect all values from this vector to indices given by localInd and store it into the output buffer...
Definition: VectorBackend.hpp:179
A container-like data-structure not storing anything and with empty implementations in many container...
Definition: FakeContainer.hpp:34
void info(int noInfoLevel, std::string const &str, Args &&... args)
prints a message, if Environment::infoLevel() >= noInfoLevel
Definition: Output.hpp:105
PetscVector(Basis const &basis)
Constructor. Constructs new BaseVector.
Definition: VectorBackend.hpp:39
PetscInt size_type
The index/size - type.
Definition: VectorBackend.hpp:34
void synchronize()
Update the ghost regions with values from the owning process.
Definition: VectorBackend.hpp:149
::Vec BaseVector
The type of the base vector.
Definition: VectorBackend.hpp:28
void finish()
Finish assembly. Must be called before vector can be used in a KSP.
Definition: VectorBackend.hpp:139
PetscScalar at(MultiIndex const &idx) const
Access the entry i of the vector with read-only-access.
Definition: VectorBackend.hpp:160
The basic container that stores a base vector data.
Definition: VectorBackend.hpp:19
static Dune::MPIHelper::MPICommunicator comm()
Return the MPI_Comm object (or a fake communicator)
Definition: Environment.hpp:86