1#ifndef DUNE_FEM_PETSCVECTOR_HH
2#define DUNE_FEM_PETSCVECTOR_HH
34 template<
class DFSpace >
53 template <
class DiscreteFunctionSpace,
class Mapper >
54 class PetscManagedDofStorage
55 :
public ManagedDofStorageImplementation< typename DiscreteFunctionSpace::GridPartType::GridType, Mapper, PetscVector< DiscreteFunctionSpace > >
57 typedef typename DiscreteFunctionSpace::GridPartType::GridType GridType;
58 typedef Mapper MapperType ;
59 typedef PetscVector< DiscreteFunctionSpace > DofArrayType ;
60 typedef ManagedDofStorageImplementation< GridType, MapperType, DofArrayType > BaseType;
65 const MapperType& mapper )
66 : BaseType( space.grid(), mapper, myArray_ ),
71 void resize(
const bool )
79 void dofCompress(
const bool clearResizedArrays )
82 if( clearResizedArrays )
89 void enableDofCompression()
91 std::cerr <<
"WARNING: PetscVector cannot handle dof compression!" << std::endl;
94 DofArrayType myArray_;
112 template<
class DFSpace >
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 >;
121 typedef PetscScalar value_type ;
122 typedef Vec DofContainerType;
124 static constexpr int blockSize = DFSpace::localBlockSize;
125 typedef Hybrid::IndexRange< int, blockSize > BlockIndices;
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;
135 typedef DofIteratorType IteratorType;
136 typedef ConstDofIteratorType ConstIteratorType;
137 typedef typename DFSpace::RangeFieldType FieldType;
138 typedef int SizeType;
140 typedef PetscMappers< DFSpace > MappersType;
142 typedef typename DFSpace::template CommDataHandle<void>::OperationType DefaultCommunicationOperationType;
145 PetscVector (
const DFSpace& space, Vec vec )
146 : mappers_( space ), vec_(vec), owner_(false)
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_ );
153 PetscVector (
const DFSpace& space )
154 : mappers_( space ), owner_(true)
156 static_assert( DefaultCommunicationOperationType::value == DFCommunicationOperation::copy ||
157 DefaultCommunicationOperationType::value == DFCommunicationOperation::add,
158 "only copy/add are available communication operations for petsc");
164 PetscVector (
const ThisType &other )
165 : mappers_( other.mappers_ ), owner_(true)
178 std::size_t size ()
const {
return mappers().ghostMapper().size(); }
180 void resize(
const std::size_t newsize = 0 )
192 void reserve(
const std::size_t capacity )
197 void hasBeenModified () { ++sequence_; }
201 communicateFlag_ =
true;
207 communicateIfNecessary();
211 const Vec* getVector ()
const
213 communicateIfNecessary();
220 communicateIfNecessary();
224 const Vec& array ()
const
226 communicateIfNecessary();
230 Vec* getGhostedVector ()
232 communicateIfNecessary();
236 const Vec* getGhostedVector ()
const
238 communicateIfNecessary();
244 ::Dune::Petsc::VecAssemblyBegin( vec_ );
249 ::Dune::Petsc::VecAssemblyEnd( vec_ );
253 template <
class Operation>
254 void communicateNow (
const Operation& operation)
const
256 communicateFlag_ =
true;
258 communicateIfNecessary( operation );
261 DofBlockType operator[] (
const IndexType index ) {
return DofBlockType( *
this,index ); }
262 ConstDofBlockType operator[] (
const IndexType index )
const {
return ConstDofBlockType( *
this,index ); }
264 ConstDofBlockPtrType block ( IndexType index )
const {
return blockPtr( index ); }
265 DofBlockPtrType block ( IndexType index ) {
return blockPtr( index ); }
267 DofBlockPtrType blockPtr ( IndexType index )
269 assert(
static_cast< std::size_t
>( index ) < mappers().size() );
270 return DofBlockPtrType(
typename DofBlockType::UnaryConstructorParamType( *
this, index ) );
273 ConstDofBlockPtrType blockPtr ( IndexType index )
const
275 assert(
static_cast< std::size_t
>( index ) < mappers().size() );
276 return ConstDofBlockPtrType(
typename ConstDofBlockType::UnaryConstructorParamType( *
this, index ) );
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() ); }
284 DofIteratorType dbegin () {
return begin(); }
285 ConstDofIteratorType dbegin ()
const {
return begin(); }
286 DofIteratorType dend () {
return end(); }
287 ConstDofIteratorType dend ()
const {
return end(); }
291 ::Dune::Petsc::VecSet( *getVector(), 0. );
292 updateGhostRegions();
293 vectorIsUpToDateNow();
296 PetscScalar operator* (
const ThisType &other )
const
299 ::Dune::Petsc::VecDot( *getVector(), *other.getVector(), &ret );
303 const ThisType& operator+= (
const ThisType &other )
305 ::Dune::Petsc::VecAXPY( *getVector(), 1., *other.getVector() );
306 updateGhostRegions();
307 vectorIsUpToDateNow();
311 const ThisType& operator-= (
const ThisType &other )
313 ::Dune::Petsc::VecAXPY( *getVector(), -1., *other.getVector() );
314 updateGhostRegions();
315 vectorIsUpToDateNow();
319 const ThisType& operator*= ( PetscScalar scalar )
321 ::Dune::Petsc::VecScale( *getVector(), scalar );
322 updateGhostRegions();
323 vectorIsUpToDateNow();
327 const ThisType& operator/= ( PetscScalar scalar )
329 assert( scalar != 0 );
330 return this->operator*=( 1./scalar );
333 void axpy (
const PetscScalar &scalar,
const ThisType &other )
335 ::Dune::Petsc::VecAXPY( *getVector(), scalar, *other.getVector() );
342 VecGetArray( ghostedVec_,&array );
343 std::fill_n( array + mappers().ghostMapper().interiorSize() * blockSize, mappers().ghostMapper().ghostSize() * blockSize, PetscScalar( 0 ) );
348 void printGlobal (
bool doit )
352 VecView( vec_, PETSC_VIEWER_STDOUT_WORLD );
355 void printGhost (
bool doit)
361 VecGetArray( ghostedVec_,&array );
362 for( std::size_t i = 0; i < size(); i++ )
364 PetscSynchronizedPrintf(PETSC_COMM_WORLD,
"%D %G\n",i,PetscRealPart(array[i]));
366 VecRestoreArray( ghostedVec_, &array );
367#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
368 PetscSynchronizedFlush( PETSC_COMM_WORLD );
370 PetscSynchronizedFlush( PETSC_COMM_WORLD, PETSC_STDOUT );
375 void assign(
const ThisType& other )
379 other.communicateIfNecessary();
382 ::Dune::Petsc::VecDuplicate( other.vec_, &vec_ );
383 ::Dune::Petsc::VecCopy( other.vec_, vec_ );
384 ::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
386 updateGhostRegions();
390 template <
class Container>
391 void assignVector(
const SimpleBlockVector< Container, blockSize >& other )
393 Vec& vec = *getGhostedVector();
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)
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 );
405 ::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
407 updateGhostRegions();
411 template <
class DofBlock>
412 void assignVector(
const ISTLBlockVector< DofBlock >& other )
414 assert( DofBlock :: dimension == blockSize );
415 Vec& vec = *getGhostedVector();
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 )
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 );
427 ::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
429 updateGhostRegions();
433 template <
class Container>
434 void copyTo( SimpleBlockVector< Container, blockSize >& other )
const
436 typedef typename Container::FieldType FieldType;
437 const PetscScalar *array =
nullptr;
438 VecGetArrayRead( ghostedVec_, &array );
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 )
446 const PetscScalar* petscBlock = array + (blockSize * mapping[ b ]);
447 FieldType* otherBlock = other.data() + (b * blockSize);
448 for(
int i=0; i<blockSize; ++i )
450 otherBlock[ i ] = petscBlock[ i ];
453 VecRestoreArrayRead( ghostedVec_, &array );
457 template <
class DofBlock>
458 void copyTo ( ISTLBlockVector< DofBlock >& other )
const
460 assert( DofBlock :: dimension == blockSize );
461 const PetscScalar *array =
nullptr;
462 VecGetArrayRead( ghostedVec_, &array );
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 )
470 const PetscScalar* petscBlock = array + (blockSize * mapping[ b ]);
471 DofBlock& otherBlock = other[ b ];
472 for(
int i=0; i<blockSize; ++i )
474 otherBlock[ i ] = petscBlock[ i ];
477 VecRestoreArrayRead( ghostedVec_, &array );
480 PetscVector& operator= (
const ThisType& other )
486 const MappersType &mappers()
const {
return mappers_; }
488 std::size_t usedMemorySize()
const
490 return size() *
sizeof(value_type);
493 static void setMemoryFactor(
const double memFactor)
499 void memMoveBackward(
const int length,
const int oldStartIdx,
const int newStartIdx)
501 DUNE_THROW(NotImplemented,
"memMoveBackward is to be implemented");
505 void memMoveForward(
const int length,
const int oldStartIdx,
const int newStartIdx)
507 DUNE_THROW(NotImplemented,
"memMoveForward is to be implemented");
510 void copyContent(
const int newIndex,
const int oldIndex )
512 DUNE_THROW(NotImplemented,
"copyContent is to be implemented");
521 const PetscInt localBlocks = mappers_.ghostMapper().interiorSize();
522 const PetscInt numGhostBlocks = mappers_.ghostMapper().ghostSize();
524 const PetscInt localSize = localBlocks * blockSize;
525 const PetscInt globalSize = mappers_.parallelMapper().size() * blockSize;
528 const PetscInt *ghostBlocks = mappers_.parallelMapper().mapping().data() + localBlocks;
529 const auto& comm = mappers_.space().gridPart().comm();
531 ::Dune::Petsc::VecCreateGhost( comm, localSize, globalSize, numGhostBlocks, ghostBlocks, &vec_ );
533 ::Dune::Petsc::VecCreateGhostBlock( comm, blockSize, localSize, globalSize, numGhostBlocks, ghostBlocks, &vec_ );
534 ::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
540 ::Dune::Petsc::VecGhostRestoreLocalForm( vec_, &ghostedVec_ );
541 ::Dune::Petsc::VecDestroy( &vec_ );
545 void communicateIfNecessary ()
const
547 DefaultCommunicationOperationType op;
548 communicateIfNecessary( op );
551 template <
class Operation>
552 void communicateIfNecessary (
const Operation& op)
const
555 if( communicateFlag_ && memorySequence_ < sequence_ )
557 if ( memorySequence_ < sequence_ )
559 if ( Operation::value == DFCommunicationOperation::add )
561 ::Dune::Petsc::VecGhostUpdateBegin( vec_, ADD_VALUES, SCATTER_REVERSE );
562 ::Dune::Petsc::VecGhostUpdateEnd( vec_, ADD_VALUES, SCATTER_REVERSE );
564 ::Dune::Petsc::VecGhostUpdateBegin( vec_, INSERT_VALUES, SCATTER_FORWARD );
565 ::Dune::Petsc::VecGhostUpdateEnd( vec_, INSERT_VALUES, SCATTER_FORWARD );
567 memorySequence_ = sequence_;
570 communicateFlag_ =
false;
575 void updateGhostRegions ()
577 ::Dune::Petsc::VecGhostUpdateBegin( vec_, INSERT_VALUES, SCATTER_FORWARD );
578 ::Dune::Petsc::VecGhostUpdateEnd( vec_, INSERT_VALUES, SCATTER_FORWARD );
581 void vectorIsUpToDateNow ()
const
583 memorySequence_ = sequence_;
584 communicateFlag_ =
false;
590 MappersType mappers_;
594 mutable unsigned long memorySequence_ = 0;
595 mutable unsigned long sequence_ = 0;
597 mutable bool communicateFlag_ =
false;
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