OGS
PETScVector.cpp
Go to the documentation of this file.
1
21#include "PETScVector.h"
22
23#include <algorithm>
24#include <cassert>
25
26#include "BaseLib/Error.h"
27
28namespace MathLib
29{
30PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size)
31{
32 if (is_global_size)
33 {
34 VecCreate(PETSC_COMM_WORLD, &v_);
35 VecSetSizes(v_, PETSC_DECIDE, vec_size);
36 }
37 else
38 {
39 // Fix size partitioning
40 // the size can be associated to specific memory allocation of a matrix
41 VecCreateMPI(PETSC_COMM_WORLD, vec_size, PETSC_DECIDE, &v_);
42 }
43
44 config();
45}
46
47PETScVector::PETScVector(const PetscInt vec_size,
48 const std::vector<PetscInt>& ghost_ids,
49 const bool is_global_size)
50 : size_ghosts_{static_cast<PetscInt>(ghost_ids.size())}, has_ghost_id_{true}
51{
52 PetscInt nghosts = static_cast<PetscInt>(ghost_ids.size());
53 if (is_global_size)
54 {
55 VecCreateGhost(PETSC_COMM_WORLD, PETSC_DECIDE, vec_size, nghosts,
56 ghost_ids.data(), &v_);
57 }
58 else
59 {
60 VecCreate(PETSC_COMM_WORLD, &v_);
61 VecSetType(v_, VECMPI);
62 VecSetSizes(v_, vec_size, PETSC_DECIDE);
63 VecMPISetGhost(v_, nghosts, ghost_ids.data());
64 }
65
66 config();
67
68 for (PetscInt i = 0; i < nghosts; i++)
69 {
70 global_ids2local_ids_ghost_.emplace(ghost_ids[i], size_loc_ + i);
71 }
72}
73
74PETScVector::PETScVector(const PETScVector& existing_vec, const bool deep_copy)
75{
76 shallowCopy(existing_vec);
77
78 // Copy values
79 if (deep_copy)
80 {
81 VecCopy(existing_vec.v_, v_);
82 }
83}
84
86 : v_{std::move(other.v_)},
87 v_loc_{std::move(other.v_loc_)},
88 start_rank_{other.start_rank_},
89 end_rank_{other.end_rank_},
90 size_{other.size_},
91 size_loc_{other.size_loc_},
92 size_ghosts_{other.size_ghosts_},
93 has_ghost_id_{other.has_ghost_id_},
94 global_ids2local_ids_ghost_{other.global_ids2local_ids_ghost_}
95{
96}
97
99{
100 VecSetFromOptions(v_);
101 VecGetOwnershipRange(v_, &start_rank_, &end_rank_);
102
103 VecGetLocalSize(v_, &size_loc_);
104 VecGetSize(v_, &size_);
105
106 VecSetOption(v_, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
107}
108
110{
111 VecAssemblyBegin(v_);
112 VecAssemblyEnd(v_);
113}
114
115void PETScVector::gatherLocalVectors(PetscScalar local_array[],
116 PetscScalar global_array[]) const
117{
118 // Collect vectors from processors.
119 int size_rank;
120 MPI_Comm_size(PETSC_COMM_WORLD, &size_rank);
121
122 // number of elements to be sent for each rank
123 std::vector<PetscInt> i_cnt(size_rank);
124
125 MPI_Allgather(&size_loc_, 1, MPI_INT, &i_cnt[0], 1, MPI_INT,
126 PETSC_COMM_WORLD);
127
128 // collect local array
129 PetscInt offset = 0;
130 // offset in the receive vector of the data from each rank
131 std::vector<PetscInt> i_disp(size_rank);
132 for (PetscInt i = 0; i < size_rank; i++)
133 {
134 i_disp[i] = offset;
135 offset += i_cnt[i];
136 }
137
138 MPI_Allgatherv(local_array, size_loc_, MPI_DOUBLE, global_array, &i_cnt[0],
139 &i_disp[0], MPI_DOUBLE, PETSC_COMM_WORLD);
140}
141
142void PETScVector::getGlobalVector(std::vector<PetscScalar>& u) const
143{
144#ifdef TEST_MEM_PETSC
145 PetscLogDouble mem1, mem2;
146 PetscMemoryGetCurrentUsage(&mem1);
147#endif
148
149 assert(static_cast<PetscInt>(u.size()) == size_);
150
151 PetscInt state;
152 VecLockGet(v_, &state);
153 if (state != 0)
154 {
155 OGS_FATAL("PETSc vector is already locked for {:s} access.",
156 state > 0 ? "read" : "write");
157 }
158 PetscScalar* xp = nullptr;
159 VecGetArray(v_, &xp);
160
161 gatherLocalVectors(xp, u.data());
162
163 // This following line may be needed late on
164 // for a communication load balance:
165 // MPI_Barrier(PETSC_COMM_WORLD);
166
167 VecRestoreArray(v_, &xp);
168
169// TEST
170#ifdef TEST_MEM_PETSC
171 PetscMemoryGetCurrentUsage(&mem2);
172 PetscPrintf(
173 PETSC_COMM_WORLD,
174 "### Memory usage by Updating. Before :%f After:%f Increase:%d\n", mem1,
175 mem2, (int)(mem2 - mem1));
176#endif
177}
178
180{
182}
183
184void PETScVector::copyValues(std::vector<PetscScalar>& u) const
185{
186 u.resize(getLocalSize() + getGhostSize());
187
188 PetscScalar* loc_x = getLocalVector();
189 std::copy_n(loc_x, getLocalSize() + getGhostSize(), u.begin());
190 restoreArray(loc_x);
191}
192
193PetscScalar PETScVector::get(const PetscInt idx) const
194{
195 if (!global_ids2local_ids_ghost_.empty())
196 {
197 return entry_array_[getLocalIndex(idx)];
198 }
199
200 // Ghost entries, and its original index is 0.
201 const PetscInt id_p = (idx == -size_) ? 0 : std::abs(idx);
202 return entry_array_[id_p];
203}
204
205std::vector<PetscScalar> PETScVector::get(
206 std::vector<IndexType> const& indices) const
207{
208 std::vector<PetscScalar> local_x(indices.size());
209 // If VecGetValues can get values from different processors,
210 // use VecGetValues(v_, indices.size(), indices.data(),
211 // local_x.data());
212
213 if (!global_ids2local_ids_ghost_.empty())
214 {
215 for (std::size_t i = 0; i < indices.size(); i++)
216 {
217 local_x[i] = entry_array_[getLocalIndex(indices[i])];
218 }
219 }
220 else
221 {
222 for (std::size_t i = 0; i < indices.size(); i++)
223 {
224 // Ghost entries, and its original index is 0.
225 const IndexType id_p =
226 (indices[i] == -size_) ? 0 : std::abs(indices[i]);
227 local_x[i] = entry_array_[id_p];
228 }
229 }
230
231 return local_x;
232}
233
234PetscScalar* PETScVector::getLocalVector() const
235{
236 PetscScalar* loc_array;
237 if (has_ghost_id_)
238 {
239 VecGhostUpdateBegin(v_, INSERT_VALUES, SCATTER_FORWARD);
240 VecGhostUpdateEnd(v_, INSERT_VALUES, SCATTER_FORWARD);
241 VecGhostGetLocalForm(v_, &v_loc_);
242 VecGetArray(v_loc_, &loc_array);
243 }
244 else
245 {
246 VecGetArray(v_, &loc_array);
247 }
248
249 return loc_array;
250}
251
252PetscInt PETScVector::getLocalIndex(const PetscInt global_index) const
253{
254 if (global_index >= 0) // non-ghost entry.
255 {
256 return global_index - start_rank_;
257 }
258
259 // A special case for a ghost location with global index equal to
260 // the size of the local vector:
261 PetscInt real_global_index = (-global_index == size_) ? 0 : -global_index;
262
263 return global_ids2local_ids_ghost_.at(real_global_index);
264}
265
266void PETScVector::restoreArray(PetscScalar* array) const
267{
268 if (has_ghost_id_)
269 {
270 VecRestoreArray(v_loc_, &array);
271 VecGhostRestoreLocalForm(v_, &v_loc_);
272 }
273 else
274 {
275 VecRestoreArray(v_, &array);
276 }
277}
278
279void PETScVector::viewer(const std::string& file_name,
280 const PetscViewerFormat vw_format) const
281{
282 PetscViewer viewer;
283 PetscViewerASCIIOpen(PETSC_COMM_WORLD, file_name.c_str(), &viewer);
284 PetscViewerPushFormat(viewer, vw_format);
285
286 PetscObjectSetName((PetscObject)v_, file_name.c_str());
287 VecView(v_, viewer);
288
289#define nEXIT_TEST
290#ifdef EXIT_TEST
291 VecDestroy(v_);
292 PetscFinalize();
293 exit(0);
294#endif
295}
296
298{
299 destroy();
300
301 VecDuplicate(v.getRawVector(), &v_);
302
303 start_rank_ = v.start_rank_;
304 end_rank_ = v.end_rank_;
305 size_ = v.size_;
306 size_loc_ = v.size_loc_;
307 size_ghosts_ = v.size_ghosts_;
308 has_ghost_id_ = v.has_ghost_id_;
309 global_ids2local_ids_ghost_ = v.global_ids2local_ids_ghost_;
310
311 config();
312}
313
315{
316 vec.finalizeAssembly();
317}
318
319} // namespace MathLib
#define OGS_FATAL(...)
Definition: Error.h:26
Declaration of class PETScVector, which provides an interface to PETSc vector routines.
Wrapper class for PETSc vector.
Definition: PETScVector.h:33
PetscInt start_rank_
Starting index in a rank.
Definition: PETScVector.h:237
void finalizeAssembly()
Perform MPI collection of assembled entries in buffer.
void gatherLocalVectors(PetscScalar local_array[], PetscScalar global_array[]) const
Collect local vectors.
std::vector< PetscScalar > entry_array_
Array containing the entries of the vector. If the vector is created without given ghost IDs,...
Definition: PETScVector.h:258
void config()
A function called by constructors to configure members.
Definition: PETScVector.cpp:98
void setLocalAccessibleVector() const
void restoreArray(PetscScalar *array) const
PetscInt getLocalSize() const
Get the number of entries in the same rank.
Definition: PETScVector.h:84
bool has_ghost_id_
Flag to indicate whether the vector is created with ghost entry indices.
Definition: PETScVector.h:249
PetscInt size_loc_
Size of local entries.
Definition: PETScVector.h:244
PetscScalar * getLocalVector() const
void getGlobalVector(std::vector< PetscScalar > &u) const
std::vector< PetscScalar > get(std::vector< IndexType > const &indices) const
PetscInt size_ghosts_
Size of local ghost entries.
Definition: PETScVector.h:246
void copyValues(std::vector< PetscScalar > &u) const
PetscInt size_
Size of the vector.
Definition: PETScVector.h:242
void shallowCopy(const PETScVector &v)
PetscInt getGhostSize() const
Get the number of ghost entries in the same rank.
Definition: PETScVector.h:86
std::map< PetscInt, PetscInt > global_ids2local_ids_ghost_
Map global indices of ghost entries to local indices.
Definition: PETScVector.h:261
PetscInt getLocalIndex(const PetscInt global_index) const
Get local index by a global index.
PetscInt end_rank_
Ending index in a rank.
Definition: PETScVector.h:239
void viewer(const std::string &file_name, const PetscViewerFormat vw_format=PETSC_VIEWER_ASCII_MATLAB) const
void finalizeVectorAssembly(VEC_T &)
General function to finalize the vector assembly.
static const double u
static const double v