AMDiS 2.10
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/common/multiindex.hh>
9
10#include <amdis/Output.hpp>
11#include <amdis/algorithm/ForEach.hpp>
12#include <amdis/algorithm/InnerProduct.hpp>
13#include <amdis/algorithm/Transform.hpp>
14#include <amdis/common/FakeContainer.hpp>
15#include <amdis/operations/Assigner.hpp>
16#include <amdis/typetree/MultiIndex.hpp>
17
18namespace AMDiS
19{
21 template <class DofMap>
23 {
24 template <class A>
25 using AssignMode = std::conditional_t<std::is_same_v<A,Assigner::plus_assign>,
26 std::integral_constant<::InsertMode, ADD_VALUES>,
27 std::integral_constant<::InsertMode, INSERT_VALUES>>;
28
29 public:
31 using BaseVector = ::Vec;
32
34 using value_type = PetscScalar;
35
37 using size_type = PetscInt;
38
39 public:
41 explicit PetscVector(DofMap const& dofMap)
42 : dofMap_(dofMap)
43 {}
44
47 : dofMap_(other.dofMap_)
48 {
49 swap(*this, other);
50 }
51
54 : dofMap_(other.dofMap_)
55 , initialized_(other.initialized_)
56 {
57 if (other.initialized_) {
58 VecDuplicate(other.vector_, &vector_);
59
60 // copy the data including the ghost values
61 Vec local, otherLocal;
62 VecGhostGetLocalForm(vector_, &local);
63 VecGhostGetLocalForm(other.vector_, &otherLocal);
64 VecCopy(otherLocal, local);
65 VecGhostRestoreLocalForm(vector_, &local);
66 VecGhostRestoreLocalForm(other.vector_, &otherLocal);
67 }
68 }
69
72 {
73 destroy();
74 }
75
78 {
79 swap(*this, other);
80 return *this;
81 }
82
85 {
86 return *this = PetscVector(other);
87 }
88
90 BaseVector const& vector() const
91 {
92 return vector_;
93 }
94
97 {
98 return vector_;
99 }
100
102 std::size_t localSize() const
103 {
104 return dofMap_.localSize();
105 }
106
108 std::size_t globalSize() const
109 {
110 return dofMap_.globalSize();
111 }
112
114 // [[possibly collective]]
115 template <class Basis>
116 void init(Basis const&, bool clear)
117 {
118 Dune::Timer t;
119 if (initialized_)
120 destroy();
121
122 // create the vector
123 create();
124 info(3, " create new vector needed {} seconds", t.elapsed());
125 t.reset();
126
127 if (clear)
128 setZero();
129
130 initialized_ = true;
131 }
132
134 // [[collective]]
135 void finish()
136 {
137 Dune::Timer t;
138 VecAssemblyBegin(vector_);
139 VecAssemblyEnd(vector_);
140 info(3, " finish vector assembly needed {} seconds", t.elapsed());
141 }
142
144 // [[collective]]
146 {
147 Dune::Timer t;
148 VecGhostUpdateBegin(vector_, INSERT_VALUES, SCATTER_FORWARD);
149 VecGhostUpdateEnd(vector_, INSERT_VALUES, SCATTER_FORWARD);
150 info(3, " synchronize vector needed {} seconds", t.elapsed());
151 }
152
154 // [[not collective]]
155 template <class MultiIndex>
156 PetscScalar at(MultiIndex const& idx) const
157 {
158 PetscScalar y = 0;
159 gather(std::array<MultiIndex,1>{idx}, &y);
160 return y;
161 }
162
164 // [[not collective]]
165 template <class MultiIndex, class Assign>
166 void insert(MultiIndex const& dof, PetscScalar value, Assign)
167 {
168 assert(initialized_);
169 VecSetValue(vector_, dofMap_.global(flatMultiIndex(dof)), value, AssignMode<Assign>::value);
170 }
171
173 // [[not collective]]
174 template <class IndexRange, class OutputIterator>
175 void gather(IndexRange const& localInd, OutputIterator buffer) const
176 {
177 forEach(localInd, [&buffer](auto&&, auto const& value) { *buffer = value; ++buffer; });
178 }
179
181
185 // [[not collective]]
186 template <class IndexRange, class LocalValues, class MaskRange, class Assign>
187 void scatter(IndexRange const& localInd, LocalValues const& localVal, MaskRange const& mask, Assign)
188 {
189 if constexpr(not std::is_same_v<typename LocalValues::value_type, PetscScalar>)
190 {
191 std::vector<PetscScalar> newLocalVal(localVal.begin(), localVal.end());
192 scatter(localInd, newLocalVal, mask, Assign{});
193 return;
194 }
195
196 assert(initialized_);
197 assert(localInd.size() == localVal.size());
198
199 // map local to global indices
200 std::vector<PetscInt> globalInd;
201 globalInd.reserve(localInd.size());
202
203 auto ind_it = std::begin(localInd);
204 auto mask_it = std::begin(mask);
205 for (; ind_it != std::end(localInd); ++mask_it, ++ind_it)
206 globalInd.push_back(*mask_it ? dofMap_.global(flatMultiIndex(*ind_it)) : -1);
207
208 if (!std::is_same_v<MaskRange, FakeContainer<bool,true>>)
209 VecSetOption(vector_, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
210
211 VecSetValues(vector_, globalInd.size(), globalInd.data(), localVal.data(), AssignMode<Assign>::value);
212 VecSetOption(vector_, VEC_IGNORE_NEGATIVE_INDICES, PETSC_FALSE);
213 }
214
216 // [[not collective]]
217 template <class IndexRange, class Func>
218 void forEach(IndexRange const& localInd, Func&& f) const
219 {
220 // Obtains the local ghosted representation of a parallel vector
221 Vec localGhostedVec = nullptr;
222 VecGhostGetLocalForm(vector_, &localGhostedVec);
223 assert(localGhostedVec != nullptr);
224
225 // A pointer to a contiguous array that contains this processor's portion of the vector data.
226 PetscScalar* ptr;
227 VecGetArray(localGhostedVec, &ptr);
228
229 for (auto const& dof : localInd)
230 {
231 size_type i = flatMultiIndex(dof);
232 if (dofMap_.owner(i)) {
233 // the index is a purely local entry
234 f(dof, ptr[dofMap_.dofToLocal(i)]);
235 }
236 else {
237 // ghost entries are stored at the end of the array
238 f(dof, ptr[dofMap_.localSize() + dofMap_.dofToGhost(i)]);
239 }
240 }
241
242 VecRestoreArray(localGhostedVec, &ptr);
243 VecGhostRestoreLocalForm(vector_, &localGhostedVec);
244 }
245
247 // [[not collective]]
248 template <class Func>
249 void forEachLocal(Func&& f) const
250 {
251 // A pointer to a contiguous array that contains this processor's portion of the vector data.
252 PetscScalar const* ptr;
253 VecGetArrayRead(vector_, &ptr);
254
255 for (size_type i = 0; i < size_type(localSize()); ++i)
256 f(ptr[dofMap_.dofToLocal(i)]);
257
258 VecRestoreArrayRead(vector_, &ptr);
259 }
260
262 // [[mutable]]
263 // [[not collective]]
264 template <class Func>
265 void forEachLocal(Func&& f)
266 {
267 // A pointer to a contiguous array that contains this processor's portion of the vector data.
268 PetscScalar *ptr;
269 VecGetArray(vector_, &ptr);
270
271 for (size_type i = 0; i < size_type(localSize()); ++i)
272 f(ptr[dofMap_.dofToLocal(i)]);
273
274 VecRestoreArray(vector_, &ptr);
275 }
276
278 // [[not collective]]
279 template <class Operation, class... VecIns>
280 static void transformLocal(PetscVector& vecOut, Operation op, VecIns const&... vecIn)
281 {
282 PetscScalar *ptrOut;
283 std::array<PetscScalar const*,sizeof...(VecIns)> ptrsIn{};
284
285 // Obtain ptrs to internal local arrays
286 VecGetArray(vecOut.vector_, &ptrOut);
287 Ranges::forIndices<sizeof...(VecIns)>([&](auto ii) {
288 VecGetArrayRead(std::get<ii>(std::tie(vecIn...)).vector_, &ptrsIn[ii]);
289 });
290
291 for (size_type i = 0; i < size_type(vecOut.localSize()); ++i) {
292 auto idx = vecOut.dofMap_.dofToLocal(i);
293 Ranges::applyIndices<sizeof...(VecIns)>([&](auto... ii) {
294 ptrOut[idx] = op(ptrsIn[ii][idx]...);
295 });
296 }
297
298 // Restore array obtained with VecGetArray[Read]()
299 VecRestoreArray(vecOut.vector_, &ptrOut);
300 Ranges::forIndices<sizeof...(VecIns)>([&](auto ii) {
301 VecRestoreArrayRead(std::get<ii>(std::tie(vecIn...)).vector_, &ptrsIn[ii]);
302 });
303 }
304
306 // [[not collective]]
307 template <class T, class BinOp1, class BinOp2>
308 static T innerProductLocal(PetscVector const& in1, PetscVector const& in2, T init, BinOp1 op1, BinOp2 op2)
309 {
310 PetscScalar const* ptr1;
311 PetscScalar const* ptr2;
312
313 // Obtain ptrs to internal local arrays
314 VecGetArrayRead(in1.vector_, &ptr1);
315 VecGetArrayRead(in2.vector_, &ptr2);
316
317 for (size_type i = 0; i < size_type(in1.localSize()); ++i) {
318 auto idx = in1.dofMap_.dofToLocal(i);
319 init = op1(std::move(init), op2(ptr1[idx], ptr2[idx]));
320 }
321
322 // Restore array obtained with VecGetArrayRead()
323 VecRestoreArrayRead(in1.vector_, &ptr1);
324 VecRestoreArrayRead(in2.vector_, &ptr2);
325
326 return init;
327 }
328
330 // [[not collective]]
331 void set(PetscScalar value)
332 {
333 Vec localGhostedVec = nullptr;
334 VecGhostGetLocalForm(vector_, &localGhostedVec);
335 assert(localGhostedVec != nullptr);
336
337 VecSet(localGhostedVec, value);
338 VecGhostRestoreLocalForm(vector_, &localGhostedVec);
339 }
340
342 // [[not collective]]
343 void setZero()
344 {
345 Vec localGhostedVec = nullptr;
346 VecGhostGetLocalForm(vector_, &localGhostedVec);
347 assert(localGhostedVec != nullptr);
348
349 VecZeroEntries(localGhostedVec);
350 VecGhostRestoreLocalForm(vector_, &localGhostedVec);
351 }
352
353 friend void swap(PetscVector& lhs, PetscVector& rhs)
354 {
355 assert(&lhs.dofMap_ == &rhs.dofMap_);
356 std::swap(lhs.vector_, rhs.vector_);
357 std::swap(lhs.initialized_, rhs.initialized_);
358 }
359
360 // An MPI Communicator or PETSC_COMM_SELF
361 MPI_Comm comm() const
362 {
363 return dofMap_.comm();
364 }
365
366 private:
367 // Create a new Vector
368 void create()
369 {
370 VecCreateGhost(comm(),
371 dofMap_.localSize(), dofMap_.globalSize(), dofMap_.ghostSize(), dofMap_.ghostIndices().data(),
372 &vector_);
373 }
374
375 // Destroy an old vector if created before
376 void destroy()
377 {
378 if (initialized_)
379 VecDestroy(&vector_);
380
381 initialized_ = false;
382 }
383
384 private:
385 // The local-to-global map
386 DofMap const& dofMap_;
387
389 mutable Vec vector_ = nullptr;
390 bool initialized_ = false;
391 };
392
393
394 namespace Recursive
395 {
396 template <class DofMap>
397 struct ForEach<PetscVector<DofMap>>
398 {
399 template <class Vec, class UnaryFunction>
400 static void impl(Vec&& vec, UnaryFunction f)
401 {
402 vec.forEachLocal(f);
403 }
404 };
405
406 template <class DofMap>
407 struct Transform<PetscVector<DofMap>>
408 {
409 template <class Operation, class... VecIns>
410 static void impl(PetscVector<DofMap>& vecOut, Operation op, VecIns const&... vecIn)
411 {
412 PetscVector<DofMap>::transformLocal(vecOut, op, vecIn...);
413 }
414 };
415
416 template <class DofMap>
417 struct InnerProduct<PetscVector<DofMap>>
418 {
419 template <class T, class BinOp1, class BinOp2>
420 static T impl(PetscVector<DofMap> const& in1, PetscVector<DofMap> const& in2, T init, BinOp1 op1, BinOp2 op2)
421 {
422 return PetscVector<DofMap>::innerProductLocal(in1, in2, std::move(init), op1, op2);
423 }
424 };
425
426 } // end namespace Recursive
427} // end namespace AMDiS
A container-like data-structure not storing anything and with empty implementations in many container...
Definition: FakeContainer.hpp:36
The basic container that stores a base vector data.
Definition: VectorBackend.hpp:23
PetscVector & operator=(PetscVector const &other)
Copy assignment operator.
Definition: VectorBackend.hpp:84
void insert(MultiIndex const &dof, PetscScalar value, Assign)
Access the entry i of the vector with write-access.
Definition: VectorBackend.hpp:166
void forEachLocal(Func &&f)
Perform the f(value) on the local elements of this vector.
Definition: VectorBackend.hpp:265
PetscVector & operator=(PetscVector &&other)
Move assignment operator.
Definition: VectorBackend.hpp:77
::Vec BaseVector
The type of the base vector.
Definition: VectorBackend.hpp:31
void forEachLocal(Func &&f) const
Perform the f(value) on the local elements of this vector.
Definition: VectorBackend.hpp:249
void setZero()
Zero all entries, including the ghost entries.
Definition: VectorBackend.hpp:343
void synchronize()
Update the ghost regions with values from the owning process.
Definition: VectorBackend.hpp:145
void finish()
Finish assembly. Must be called before vector can be used in a KSP.
Definition: VectorBackend.hpp:135
PetscVector(PetscVector &&other)
Move constructor.
Definition: VectorBackend.hpp:46
PetscScalar at(MultiIndex const &idx) const
Access the entry i of the vector with read-only-access.
Definition: VectorBackend.hpp:156
static T innerProductLocal(PetscVector const &in1, PetscVector const &in2, T init, BinOp1 op1, BinOp2 op2)
Perform the op1(init, op2(value1, value2)) on the local elements of [in1] and [in2].
Definition: VectorBackend.hpp:308
PetscVector(PetscVector const &other)
Copy constructor.
Definition: VectorBackend.hpp:53
~PetscVector()
Destructor.
Definition: VectorBackend.hpp:71
BaseVector & vector()
Return the data-vector vector_.
Definition: VectorBackend.hpp:96
std::size_t globalSize() const
Return the number of entries in the global vector.
Definition: VectorBackend.hpp:108
static void transformLocal(PetscVector &vecOut, Operation op, VecIns const &... vecIn)
Perform the valueOut = op(valueIn...) on the local elements of [vecOut] and [vecIn....
Definition: VectorBackend.hpp:280
void init(Basis const &, bool clear)
Resize the vector_ to the size s.
Definition: VectorBackend.hpp:116
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:187
PetscVector(DofMap const &dofMap)
Constructor. Constructs new BaseVector.
Definition: VectorBackend.hpp:41
std::size_t localSize() const
Return the number of entries in the local part of the vector.
Definition: VectorBackend.hpp:102
PetscInt size_type
The index/size - type.
Definition: VectorBackend.hpp:37
BaseVector const & vector() const
Return the data-vector vector_.
Definition: VectorBackend.hpp:90
void set(PetscScalar value)
Set all entries to value, including the ghost entries.
Definition: VectorBackend.hpp:331
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:175
PetscScalar value_type
The type of the elements of the DOFVector.
Definition: VectorBackend.hpp:34
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:218
constexpr bool MultiIndex
A multi-index type.
Definition: Concepts.hpp:149
Definition: ForEach.hpp:39
General implementation of recursive inner-product.
Definition: InnerProduct.hpp:58
Definition: Transform.hpp:35