dune-fem 2.8.0
Loading...
Searching...
No Matches
petscvector.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_PETSCVECTOR_HH
2#define DUNE_FEM_PETSCVECTOR_HH
3
4#include <cassert>
5#include <cstddef>
6
7#include <algorithm>
8#include <string>
9
12
14
16
17#if HAVE_PETSC
18
24
25namespace Dune
26{
27
28 namespace Fem
29 {
30
31 // Internal Forward Declarations
32 // -----------------------------
33
34 template< class DFSpace >
35 class PetscVector;
36
37
38 // External Forward Declarations
39 // -----------------------------
40
41 template< class >
42 class PetscDofBlock;
43
44 template< class >
45 class PetscDofProxy;
46
47
48
49 // PetscManagedDofStorage
50 // ----------------------
51
53 template < class DiscreteFunctionSpace, class Mapper >
54 class PetscManagedDofStorage
55 : public ManagedDofStorageImplementation< typename DiscreteFunctionSpace::GridPartType::GridType, Mapper, PetscVector< DiscreteFunctionSpace > >
56 {
57 typedef typename DiscreteFunctionSpace::GridPartType::GridType GridType;
58 typedef Mapper MapperType ;
59 typedef PetscVector< DiscreteFunctionSpace > DofArrayType ;
60 typedef ManagedDofStorageImplementation< GridType, MapperType, DofArrayType > BaseType;
61
62 public:
64 PetscManagedDofStorage( const DiscreteFunctionSpace& space,
65 const MapperType& mapper )
66 : BaseType( space.grid(), mapper, myArray_ ),
67 myArray_( space )
68 {
69 }
70
71 void resize( const bool )
72 { // do nothing here, only in compress
73 }
74 void reserve( int )
75 {
76 // do nothing here, only in compress
77 }
78
79 void dofCompress( const bool clearResizedArrays )
80 {
81 myArray_.resize();
82 if( clearResizedArrays )
83 {
84 myArray_.clear();
85 }
86 }
87
89 void enableDofCompression()
90 {
91 std::cerr << "WARNING: PetscVector cannot handle dof compression!" << std::endl;
92 }
93 protected:
94 DofArrayType myArray_;
95 };
96
97
98 // PetscVector
99 // -----------
100
101 /*
102 * This encapsules a PETSc Vec with ghosts.
103 * Some conceptual explanations:
104 * The PETSc vector, as modelled by this class, consists of three parts:
105 * 1) the whole PETSc vector, which might be distributed across several MPI processes.
106 * We call this the _global vector_
107 * 2) Each process holds a portion of this global vector, we call this part the
108 * _process vector_.
109 * 3) And there is a representation of the process vector, which also has 'ghost dof blocks'.
110 * We call this represenation the _ghosted vector_.
111 */
112 template< class DFSpace >
113 class PetscVector
114 {
115 typedef PetscVector< DFSpace > ThisType;
116 friend class PetscDofBlock< ThisType >;
117 friend class PetscDofBlock< const ThisType >;
118 friend class PetscDofProxy< ThisType >;
119 friend class PetscDofProxy< const ThisType >;
120 public:
121 typedef PetscScalar value_type ;
122 typedef Vec DofContainerType;
123
124 static constexpr int blockSize = DFSpace::localBlockSize;
125 typedef Hybrid::IndexRange< int, blockSize > BlockIndices;
126
127 typedef PetscDofBlock< ThisType > DofBlockType;
128 typedef PetscDofBlock< const ThisType > ConstDofBlockType;
129 typedef typename DofBlockType::DofIterator DofIteratorType;
130 typedef typename ConstDofBlockType::DofIterator ConstDofIteratorType;
131 typedef Envelope< DofBlockType > DofBlockPtrType;
132 typedef Envelope< ConstDofBlockType > ConstDofBlockPtrType;
133 typedef typename DofBlockType::IndexType IndexType;
134
135 typedef DofIteratorType IteratorType;
136 typedef ConstDofIteratorType ConstIteratorType;
137 typedef typename DFSpace::RangeFieldType FieldType;
138 typedef int SizeType;
139
140 typedef PetscMappers< DFSpace > MappersType;
141
142 typedef typename DFSpace::template CommDataHandle<void>::OperationType DefaultCommunicationOperationType;
143
144 // note that Vec is a pointer type so no deep copy is made
145 PetscVector ( const DFSpace& space, Vec vec )
146 : mappers_( space ), vec_(vec), owner_(false)
147 {
148 static_assert( DefaultCommunicationOperationType::value == DFCommunicationOperation::copy ||
149 DefaultCommunicationOperationType::value == DFCommunicationOperation::add,
150 "only copy/add are available communication operations for petsc");
151 ::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
152 }
153 PetscVector ( const DFSpace& space )
154 : mappers_( space ), owner_(true)
155 {
156 static_assert( DefaultCommunicationOperationType::value == DFCommunicationOperation::copy ||
157 DefaultCommunicationOperationType::value == DFCommunicationOperation::add,
158 "only copy/add are available communication operations for petsc");
159 // init vector
160 init();
161 }
162
163 // TODO: think about sequence overflows...
164 PetscVector ( const ThisType &other )
165 : mappers_( other.mappers_ ), owner_(true)
166 {
167 // assign vectors
168 assign( other );
169 }
170
171 ~PetscVector ()
172 {
173 if (owner_)
174 // destroy vectors
175 removeObj();
176 }
177
178 std::size_t size () const { return mappers().ghostMapper().size(); }
179
180 void resize( const std::size_t newsize = 0 )
181 {
182 // TODO: keep old data stored in current vector
183 // remove old vectors
184 removeObj();
185
186 // initialize new
187 init();
188
189 hasBeenModified ();
190 }
191
192 void reserve( const std::size_t capacity )
193 {
194 resize( capacity );
195 }
196
197 void hasBeenModified () { ++sequence_; }
198
199 void communicate ()
200 {
201 communicateFlag_ = true;
202 }
203
204 // accessors for the underlying PETSc vectors
205 Vec* getVector ()
206 {
207 communicateIfNecessary();
208 return &vec_;
209 }
210
211 const Vec* getVector () const
212 {
213 communicateIfNecessary();
214 return &vec_;
215 }
216
217 // accessors for the underlying PETSc vectors
218 Vec& array ()
219 {
220 communicateIfNecessary();
221 return vec_;
222 }
223
224 const Vec& array () const
225 {
226 communicateIfNecessary();
227 return vec_;
228 }
229
230 Vec* getGhostedVector ()
231 {
232 communicateIfNecessary();
233 return &ghostedVec_;
234 }
235
236 const Vec* getGhostedVector () const
237 {
238 communicateIfNecessary();
239 return &ghostedVec_;
240 }
241
242 void beginAssemble()
243 {
244 ::Dune::Petsc::VecAssemblyBegin( vec_ );
245 }
246
247 void endAssemble()
248 {
249 ::Dune::Petsc::VecAssemblyEnd( vec_ );
250 }
251
252 // force communication _now_
253 template <class Operation>
254 void communicateNow (const Operation& operation) const
255 {
256 communicateFlag_ = true;
257 ++sequence_;
258 communicateIfNecessary( operation );
259 }
260
261 DofBlockType operator[] ( const IndexType index ) { return DofBlockType( *this,index ); }
262 ConstDofBlockType operator[] ( const IndexType index ) const { return ConstDofBlockType( *this,index ); }
263
264 ConstDofBlockPtrType block ( IndexType index ) const { return blockPtr( index ); }
265 DofBlockPtrType block ( IndexType index ) { return blockPtr( index ); }
266
267 DofBlockPtrType blockPtr ( IndexType index )
268 {
269 assert( static_cast< std::size_t >( index ) < mappers().size() );
270 return DofBlockPtrType( typename DofBlockType::UnaryConstructorParamType( *this, index ) );
271 }
272
273 ConstDofBlockPtrType blockPtr ( IndexType index ) const
274 {
275 assert( static_cast< std::size_t >( index ) < mappers().size() );
276 return ConstDofBlockPtrType( typename ConstDofBlockType::UnaryConstructorParamType( *this, index ) );
277 }
278
279 DofIteratorType begin () { return DofIteratorType( *this, 0, 0 ); }
280 ConstDofIteratorType begin () const { return ConstDofIteratorType( *this, 0, 0 ); }
281 DofIteratorType end () { return DofIteratorType( *this, mappers().size() ); }
282 ConstDofIteratorType end () const { return ConstDofIteratorType( *this, mappers().size() ); }
283
284 DofIteratorType dbegin () { return begin(); }
285 ConstDofIteratorType dbegin () const { return begin(); }
286 DofIteratorType dend () { return end(); }
287 ConstDofIteratorType dend () const { return end(); }
288
289 void clear ()
290 {
291 ::Dune::Petsc::VecSet( *getVector(), 0. );
292 updateGhostRegions();
293 vectorIsUpToDateNow();
294 }
295
296 PetscScalar operator* ( const ThisType &other ) const
297 {
298 PetscScalar ret;
299 ::Dune::Petsc::VecDot( *getVector(), *other.getVector(), &ret );
300 return ret;
301 }
302
303 const ThisType& operator+= ( const ThisType &other )
304 {
305 ::Dune::Petsc::VecAXPY( *getVector(), 1., *other.getVector() );
306 updateGhostRegions();
307 vectorIsUpToDateNow();
308 return *this;
309 }
310
311 const ThisType& operator-= ( const ThisType &other )
312 {
313 ::Dune::Petsc::VecAXPY( *getVector(), -1., *other.getVector() );
314 updateGhostRegions();
315 vectorIsUpToDateNow();
316 return *this;
317 }
318
319 const ThisType& operator*= ( PetscScalar scalar )
320 {
321 ::Dune::Petsc::VecScale( *getVector(), scalar );
322 updateGhostRegions();
323 vectorIsUpToDateNow();
324 return *this;
325 }
326
327 const ThisType& operator/= ( PetscScalar scalar )
328 {
329 assert( scalar != 0 );
330 return this->operator*=( 1./scalar );
331 }
332
333 void axpy ( const PetscScalar &scalar, const ThisType &other )
334 {
335 ::Dune::Petsc::VecAXPY( *getVector(), scalar, *other.getVector() );
336 hasBeenModified();
337 }
338
339 void clearGhost( )
340 {
341 PetscScalar *array;
342 VecGetArray( ghostedVec_,&array );
343 std::fill_n( array + mappers().ghostMapper().interiorSize() * blockSize, mappers().ghostMapper().ghostSize() * blockSize, PetscScalar( 0 ) );
344 }
345
346 // debugging; comes in handy to call these 2 methods in gdb
347 // doit is only here to prevent the compiler from optimizing these calls away...
348 void printGlobal ( bool doit )
349 {
350 if( !doit )
351 return;
352 VecView( vec_, PETSC_VIEWER_STDOUT_WORLD );
353 }
354
355 void printGhost ( bool doit)
356 {
357 if( !doit )
358 return;
359
360 PetscScalar *array;
361 VecGetArray( ghostedVec_,&array );
362 for( std::size_t i = 0; i < size(); i++ )
363 {
364 PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%D %G\n",i,PetscRealPart(array[i]));
365 }
366 VecRestoreArray( ghostedVec_, &array );
367#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
368 PetscSynchronizedFlush( PETSC_COMM_WORLD );
369#else
370 PetscSynchronizedFlush( PETSC_COMM_WORLD, PETSC_STDOUT );
371#endif
372 }
373
374 // assign from other given PetscVector
375 void assign( const ThisType& other )
376 {
377 // we want the 'other' to do all its communication right now before
378 // we start copying values from it
379 other.communicateIfNecessary();
380
381 // Do the copying on the PETSc level
382 ::Dune::Petsc::VecDuplicate( other.vec_, &vec_ );
383 ::Dune::Petsc::VecCopy( other.vec_, vec_ );
384 ::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
385
386 updateGhostRegions();
387 }
388
389 // assign from other given SimpleBlockVector with same block size
390 template <class Container>
391 void assignVector( const SimpleBlockVector< Container, blockSize >& other )
392 {
393 Vec& vec = *getGhostedVector();
394
395 // use mapping from the ghost mapper which removes duplicate dofs
396 const auto& mapping = mappers_.ghostMapper().mapping();
397 const size_t blocks = mapping.size();
398 assert( blocks == other.size() );
399 for( size_t b=0, bs = 0; b<blocks; ++b, bs += blockSize)
400 {
401 PetscInt block = mapping[ b ];
402 const PetscScalar* values = static_cast< const PetscScalar* > (other.data()+bs);
403 ::Dune::Petsc::VecSetValuesBlocked( vec, 1, &block, values, INSERT_VALUES );
404 }
405 ::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
406
407 updateGhostRegions();
408 }
409
410 // assign from other given ISTLBlockVector with same block size
411 template <class DofBlock>
412 void assignVector( const ISTLBlockVector< DofBlock >& other )
413 {
414 assert( DofBlock :: dimension == blockSize );
415 Vec& vec = *getGhostedVector();
416
417 // use mapping from the ghost mapper which removes duplicate dofs
418 const auto& mapping = mappers_.ghostMapper().mapping();
419 const size_t blocks = mapping.size();
420 assert( blocks == other.size() );
421 for( size_t b=0; b<blocks; ++b )
422 {
423 PetscInt block = mapping[ b ];
424 const PetscScalar* values = static_cast< const PetscScalar* > (&(other[ b ][ 0 ])) ;
425 ::Dune::Petsc::VecSetValuesBlocked( vec, 1, &block, values, INSERT_VALUES );
426 }
427 ::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
428
429 updateGhostRegions();
430 }
431
432 // assign from other given SimpleBlockVector with same block size
433 template <class Container>
434 void copyTo( SimpleBlockVector< Container, blockSize >& other ) const
435 {
436 typedef typename Container::FieldType FieldType;
437 const PetscScalar *array = nullptr;
438 VecGetArrayRead( ghostedVec_, &array );
439
440 // use mapping from the ghost mapper which removes duplicate dofs
441 const auto& mapping = mappers_.ghostMapper().mapping();
442 const size_t blocks = mapping.size();
443 assert( blocks == other.size() );
444 for( size_t b=0; b<blocks; ++b )
445 {
446 const PetscScalar* petscBlock = array + (blockSize * mapping[ b ]);
447 FieldType* otherBlock = other.data() + (b * blockSize);
448 for( int i=0; i<blockSize; ++i )
449 {
450 otherBlock[ i ] = petscBlock[ i ];
451 }
452 }
453 VecRestoreArrayRead( ghostedVec_, &array );
454 }
455
456 // assign from other given ISTLBlockVector with same block size
457 template <class DofBlock>
458 void copyTo ( ISTLBlockVector< DofBlock >& other ) const
459 {
460 assert( DofBlock :: dimension == blockSize );
461 const PetscScalar *array = nullptr;
462 VecGetArrayRead( ghostedVec_, &array );
463
464 // use mapping from the ghost mapper which removes duplicate dofs
465 const auto& mapping = mappers_.ghostMapper().mapping();
466 const size_t blocks = mapping.size();
467 assert( blocks == other.size() );
468 for( size_t b=0; b<blocks; ++b )
469 {
470 const PetscScalar* petscBlock = array + (blockSize * mapping[ b ]);
471 DofBlock& otherBlock = other[ b ];
472 for( int i=0; i<blockSize; ++i )
473 {
474 otherBlock[ i ] = petscBlock[ i ];
475 }
476 }
477 VecRestoreArrayRead( ghostedVec_, &array );
478 }
479
480 PetscVector& operator= ( const ThisType& other )
481 {
482 assign( other );
483 return *this;
484 }
485
486 const MappersType &mappers() const { return mappers_; }
487
488 std::size_t usedMemorySize() const
489 {
490 return size() * sizeof(value_type);
491 }
492
493 static void setMemoryFactor(const double memFactor)
494 {
495 // do nothing here
496 }
497
499 void memMoveBackward( const int length, const int oldStartIdx, const int newStartIdx)
500 {
501 DUNE_THROW(NotImplemented,"memMoveBackward is to be implemented");
502 }
503
505 void memMoveForward(const int length, const int oldStartIdx, const int newStartIdx)
506 {
507 DUNE_THROW(NotImplemented,"memMoveForward is to be implemented");
508 }
509
510 void copyContent( const int newIndex, const int oldIndex )
511 {
512 DUNE_THROW(NotImplemented,"copyContent is to be implemented");
513 }
514
515 protected:
516 // setup vector according to mapping sizes
517 void init()
518 {
519 mappers_.update();
520
521 const PetscInt localBlocks = mappers_.ghostMapper().interiorSize();
522 const PetscInt numGhostBlocks = mappers_.ghostMapper().ghostSize();
523
524 const PetscInt localSize = localBlocks * blockSize;
525 const PetscInt globalSize = mappers_.parallelMapper().size() * blockSize;
526
527 // finally, create the PETSc Vecs
528 const PetscInt *ghostBlocks = mappers_.parallelMapper().mapping().data() + localBlocks;
529 const auto& comm = mappers_.space().gridPart().comm();
530 if( blockSize == 1 )
531 ::Dune::Petsc::VecCreateGhost( comm, localSize, globalSize, numGhostBlocks, ghostBlocks, &vec_ );
532 else
533 ::Dune::Petsc::VecCreateGhostBlock( comm, blockSize, localSize, globalSize, numGhostBlocks, ghostBlocks, &vec_ );
534 ::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
535 }
536
537 // delete vectors
538 void removeObj()
539 {
540 ::Dune::Petsc::VecGhostRestoreLocalForm( vec_, &ghostedVec_ );
541 ::Dune::Petsc::VecDestroy( &vec_ );
542 }
543
544 // communicate using the space default communication option
545 void communicateIfNecessary () const
546 {
547 DefaultCommunicationOperationType op;
548 communicateIfNecessary( op );
549 }
550
551 template <class Operation>
552 void communicateIfNecessary (const Operation& op) const
553 {
554 // communicate this process' values
555 if( communicateFlag_ && memorySequence_ < sequence_ )
556 {
557 if ( memorySequence_ < sequence_ )
558 {
559 if ( Operation::value == DFCommunicationOperation::add )
560 {
561 ::Dune::Petsc::VecGhostUpdateBegin( vec_, ADD_VALUES, SCATTER_REVERSE );
562 ::Dune::Petsc::VecGhostUpdateEnd( vec_, ADD_VALUES, SCATTER_REVERSE );
563 }
564 ::Dune::Petsc::VecGhostUpdateBegin( vec_, INSERT_VALUES, SCATTER_FORWARD );
565 ::Dune::Petsc::VecGhostUpdateEnd( vec_, INSERT_VALUES, SCATTER_FORWARD );
566
567 memorySequence_ = sequence_;
568 }
569
570 communicateFlag_ = false;
571 }
572 }
573
574 // Updates the ghost dofs, obtains them from the owning process
575 void updateGhostRegions ()
576 {
577 ::Dune::Petsc::VecGhostUpdateBegin( vec_, INSERT_VALUES, SCATTER_FORWARD );
578 ::Dune::Petsc::VecGhostUpdateEnd( vec_, INSERT_VALUES, SCATTER_FORWARD );
579 }
580
581 void vectorIsUpToDateNow () const
582 {
583 memorySequence_ = sequence_;
584 communicateFlag_ = false;
585 }
586
587 /*
588 * data fields
589 */
590 MappersType mappers_;
591 Vec vec_;
592 Vec ghostedVec_;
593
594 mutable unsigned long memorySequence_ = 0; // represents the state of the PETSc vec in memory
595 mutable unsigned long sequence_ = 0; // represents the modifications to the PETSc vec
596
597 mutable bool communicateFlag_ = false;
598 bool owner_;
599 };
600
601 } // namespace Fem
602
603} // namespace Dune
604
605#endif // #if HAVE_PETSC
606
607#endif // DUNE_FEM_PETSCVECTOR_HH
Definition: bindguard.hh:11
void axpy(const T &a, const T &x, T &y)
Definition: space/basisfunctionset/functor.hh:38
typename Impl::Index< Range >::Type IndexType
Definition: hybrid.hh:69
discrete function space