dune-fem  2.8-git
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 
25 namespace 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
@ copy
Definition: commoperations.hh:123
@ add
Definition: commoperations.hh:123
Definition: bindguard.hh:11
Double operator*(const Double &a, const Double &b)
Definition: double.hh:506
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