00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
00021 #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
00022
00023 #include <opm/common/Exceptions.hpp>
00024
00025 #include <dune/istl/preconditioner.hh>
00026 #include <dune/istl/paamg/smoother.hh>
00027 #include <dune/istl/paamg/pinfo.hh>
00028
00029 namespace Opm
00030 {
00031
00032
00033
00034 template<class Matrix, class Domain, class Range, class ParallelInfo = Dune::Amg::SequentialInformation>
00035 class ParallelOverlappingILU0;
00036
00037 }
00038
00039 namespace Dune
00040 {
00041
00042 namespace Amg
00043 {
00044
00051 template<class Matrix, class Domain, class Range, class ParallelInfo>
00052 struct ConstructionTraits<Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> >
00053 {
00054 typedef Dune::SeqILU0<Matrix,Domain,Range> T;
00055 typedef DefaultParallelConstructionArgs<T,ParallelInfo> Arguments;
00056 typedef ConstructionTraits<T> SeqConstructionTraits;
00057 static inline Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* construct(Arguments& args)
00058 {
00059 return new Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>(args.getMatrix(),
00060 args.getComm(),
00061 args.getArgs().relaxationFactor);
00062 }
00063
00064 static inline void deconstruct(Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* bp)
00065 {
00066 delete bp;
00067 }
00068
00069 };
00070
00071 }
00072
00073
00074
00075 }
00076
00077 namespace Opm
00078 {
00079 namespace detail
00080 {
00082 template<class M, class CRS, class InvVector>
00083 void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv )
00084 {
00085 typedef typename M :: size_type size_type;
00086
00087 lower.resize( A.N() );
00088 upper.resize( A.N() );
00089 inv.resize( A.N() );
00090
00091 lower.reserveAdditional( 2*A.N() );
00092
00093
00094 const auto endi = A.end();
00095 size_type row = 0;
00096 size_type colcount = 0;
00097 lower.rows_[ 0 ] = colcount;
00098 for (auto i=A.begin(); i!=endi; ++i, ++row)
00099 {
00100 const size_type iIndex = i.index();
00101 lower.reserveAdditional( (*i).size() );
00102
00103
00104 for (auto j=(*i).begin(); j.index() < iIndex; ++j )
00105 {
00106 lower.push_back( (*j), j.index() );
00107 ++colcount;
00108 }
00109 lower.rows_[ iIndex+1 ] = colcount;
00110 }
00111
00112 const auto rendi = A.beforeBegin();
00113 row = 0;
00114 colcount = 0;
00115 upper.rows_[ 0 ] = colcount ;
00116
00117 upper.reserveAdditional( lower.nonZeros() + A.N() );
00118
00119
00120
00121 for (auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
00122 {
00123 const size_type iIndex = i.index();
00124 upper.reserveAdditional( (*i).size() );
00125
00126
00127
00128 for (auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
00129 {
00130 const size_type jIndex = j.index();
00131 if( j.index() == iIndex )
00132 {
00133 inv[ row ] = (*j);
00134 break;
00135 }
00136 else if ( j.index() >= i.index() )
00137 {
00138 upper.push_back( (*j), jIndex );
00139 ++colcount ;
00140 }
00141 }
00142 upper.rows_[ row+1 ] = colcount;
00143 }
00144 }
00145 }
00146
00162 template<class Matrix, class Domain, class Range, class ParallelInfoT>
00163 class ParallelOverlappingILU0
00164 : public Dune::Preconditioner<Domain,Range>
00165 {
00166 typedef ParallelInfoT ParallelInfo;
00167
00168
00169 public:
00171 typedef typename std::remove_const<Matrix>::type matrix_type;
00173 typedef Domain domain_type;
00175 typedef Range range_type;
00177 typedef typename Domain::field_type field_type;
00178
00179 typedef typename matrix_type::block_type block_type;
00180 typedef typename matrix_type::size_type size_type;
00181
00182 protected:
00183 struct CRS
00184 {
00185 CRS() : nRows_( 0 ) {}
00186
00187 size_type rows() const { return nRows_; }
00188
00189 size_type nonZeros() const
00190 {
00191 assert( rows_[ rows() ] != size_type(-1) );
00192 return rows_[ rows() ];
00193 }
00194
00195 void resize( const size_type nRows )
00196 {
00197 if( nRows_ != nRows )
00198 {
00199 nRows_ = nRows ;
00200 rows_.resize( nRows_+1, size_type(-1) );
00201 }
00202 }
00203
00204 void reserveAdditional( const size_type nonZeros )
00205 {
00206 const size_type needed = values_.size() + nonZeros ;
00207 if( values_.capacity() < needed )
00208 {
00209 const size_type estimate = needed * 1.1;
00210 values_.reserve( estimate );
00211 cols_.reserve( estimate );
00212 }
00213 }
00214
00215 void push_back( const block_type& value, const size_type index )
00216 {
00217 values_.push_back( value );
00218 cols_.push_back( index );
00219 }
00220
00221 std::vector< size_type > rows_;
00222 std::vector< block_type > values_;
00223 std::vector< size_type > cols_;
00224 size_type nRows_;
00225 };
00226
00227 public:
00228
00229 enum {
00231 category = std::is_same<ParallelInfoT, Dune::Amg::SequentialInformation>::value ?
00232 Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping
00233 };
00234
00242 template<class BlockType, class Alloc>
00243 ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
00244 const int n, const field_type w )
00245 : lower_(),
00246 upper_(),
00247 inv_(),
00248 comm_(nullptr), w_(w),
00249 relaxation_( std::abs( w - 1.0 ) > 1e-15 )
00250 {
00251
00252
00253 init( reinterpret_cast<const Matrix&>(A), n );
00254 }
00255
00262 template<class BlockType, class Alloc>
00263 ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
00264 const field_type w)
00265 : ParallelOverlappingILU0( A, 0, w )
00266 {
00267 }
00268
00276 template<class BlockType, class Alloc>
00277 ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
00278 const ParallelInfo& comm, const field_type w)
00279 : lower_(),
00280 upper_(),
00281 inv_(),
00282 comm_(&comm), w_(w),
00283 relaxation_( std::abs( w - 1.0 ) > 1e-15 )
00284 {
00285
00286
00287 init( reinterpret_cast<const Matrix&>(A), 0 );
00288 }
00289
00295 virtual void pre (Domain& x, Range& b)
00296 {
00297 DUNE_UNUSED_PARAMETER(x);
00298 DUNE_UNUSED_PARAMETER(b);
00299 }
00300
00306 virtual void apply (Domain& v, const Range& d)
00307 {
00308 Range& md = const_cast<Range&>(d);
00309 copyOwnerToAll( md );
00310
00311
00312 typedef typename Range ::block_type dblock;
00313 typedef typename Domain::block_type vblock;
00314
00315 const size_type iEnd = lower_.rows();
00316 const size_type lastRow = iEnd - 1;
00317 if( iEnd != upper_.rows() )
00318 {
00319 std::abort();
00320
00321 }
00322
00323
00324 for( size_type i=0; i<iEnd; ++ i )
00325 {
00326 dblock rhs( d[ i ] );
00327 const size_type rowI = lower_.rows_[ i ];
00328 const size_type rowINext = lower_.rows_[ i+1 ];
00329
00330 for( size_type col = rowI; col < rowINext; ++ col )
00331 {
00332 lower_.values_[ col ].mmv( v[ lower_.cols_[ col ] ], rhs );
00333 }
00334
00335 v[ i ] = rhs;
00336 }
00337
00338 copyOwnerToAll( v );
00339
00340 for( size_type i=0; i<iEnd; ++ i )
00341 {
00342 vblock& vBlock = v[ lastRow - i ];
00343 vblock rhs ( vBlock );
00344 const size_type rowI = upper_.rows_[ i ];
00345 const size_type rowINext = upper_.rows_[ i+1 ];
00346
00347 for( size_type col = rowI; col < rowINext; ++ col )
00348 {
00349 upper_.values_[ col ].mmv( v[ upper_.cols_[ col ] ], rhs );
00350 }
00351
00352
00353 inv_[ i ].mv( rhs, vBlock);
00354 }
00355
00356 copyOwnerToAll( v );
00357
00358 if( relaxation_ ) {
00359 v *= w_;
00360 }
00361 }
00362
00363 template <class V>
00364 void copyOwnerToAll( V& v ) const
00365 {
00366 if( comm_ ) {
00367 comm_->copyOwnerToAll(v, v);
00368 }
00369 }
00370
00376 virtual void post (Range& x)
00377 {
00378 DUNE_UNUSED_PARAMETER(x);
00379 }
00380
00381 protected:
00382 void init( const Matrix& A, const int iluIteration )
00383 {
00384 int ilu_setup_successful = 1;
00385 std::string message;
00386 const int rank = ( comm_ ) ? comm_->communicator().rank() : 0;
00387
00388 std::unique_ptr< Matrix > ILU;
00389
00390 try
00391 {
00392 if( iluIteration == 0 ) {
00393
00394 ILU.reset( new Matrix( A ) );
00395 bilu0_decomposition( *ILU );
00396 }
00397 else {
00398
00399 ILU.reset( new Matrix( A.N(), A.M(), Matrix::row_wise) );
00400 bilu_decomposition( A, iluIteration, *ILU );
00401 }
00402 }
00403 catch ( Dune::MatrixBlockError error )
00404 {
00405 message = error.what();
00406 std::cerr<<"Exception occured on process " << rank << " during " <<
00407 "setup of ILU0 preconditioner with message: " <<
00408 message<<std::endl;
00409 ilu_setup_successful = 0;
00410 }
00411
00412
00413 if ( comm_ && comm_->communicator().min(ilu_setup_successful) == 0 )
00414 {
00415 throw Dune::MatrixBlockError();
00416 }
00417
00418
00419 detail::convertToCRS( *ILU, lower_, upper_, inv_ );
00420 }
00421
00422 protected:
00424 CRS lower_;
00425 CRS upper_;
00426 std::vector< block_type > inv_;
00427
00428 const ParallelInfo* comm_;
00430 const field_type w_;
00431 const bool relaxation_;
00432
00433 };
00434
00435 }
00436 #endif