20 #ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
21 #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
23 #include <opm/common/Exceptions.hpp>
25 #include <dune/istl/preconditioner.hh>
26 #include <dune/istl/paamg/smoother.hh>
27 #include <dune/istl/paamg/pinfo.hh>
34 template<
class Matrix,
class Domain,
class Range,
class ParallelInfo = Dune::Amg::SequentialInformation>
51 template<
class Matrix,
class Domain,
class Range,
class ParallelInfo>
52 struct ConstructionTraits<Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> >
54 typedef Dune::SeqILU0<Matrix,Domain,Range> T;
55 typedef DefaultParallelConstructionArgs<T,ParallelInfo> Arguments;
56 typedef ConstructionTraits<T> SeqConstructionTraits;
61 args.getArgs().relaxationFactor);
82 template<
class M,
class CRS,
class InvVector>
83 void convertToCRS(
const M& A, CRS& lower, CRS& upper, InvVector& inv )
85 typedef typename M :: size_type size_type;
87 lower.resize( A.N() );
88 upper.resize( A.N() );
91 lower.reserveAdditional( 2*A.N() );
94 const auto endi = A.end();
96 size_type colcount = 0;
97 lower.rows_[ 0 ] = colcount;
98 for (
auto i=A.begin(); i!=endi; ++i, ++row)
100 const size_type iIndex = i.index();
101 lower.reserveAdditional( (*i).size() );
104 for (
auto j=(*i).begin(); j.index() < iIndex; ++j )
106 lower.push_back( (*j), j.index() );
109 lower.rows_[ iIndex+1 ] = colcount;
112 const auto rendi = A.beforeBegin();
115 upper.rows_[ 0 ] = colcount ;
117 upper.reserveAdditional( lower.nonZeros() + A.N() );
121 for (
auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
123 const size_type iIndex = i.index();
124 upper.reserveAdditional( (*i).size() );
128 for (
auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
130 const size_type jIndex = j.index();
131 if( j.index() == iIndex )
136 else if ( j.index() >= i.index() )
138 upper.push_back( (*j), jIndex );
142 upper.rows_[ row+1 ] = colcount;
162 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
163 class ParallelOverlappingILU0
164 :
public Dune::Preconditioner<Domain,Range>
166 typedef ParallelInfoT ParallelInfo;
179 typedef typename matrix_type::block_type block_type;
180 typedef typename matrix_type::size_type size_type;
185 CRS() : nRows_( 0 ) {}
187 size_type rows()
const {
return nRows_; }
189 size_type nonZeros()
const
191 assert( rows_[ rows() ] != size_type(-1) );
192 return rows_[ rows() ];
195 void resize(
const size_type nRows )
197 if( nRows_ != nRows )
200 rows_.resize( nRows_+1, size_type(-1) );
204 void reserveAdditional(
const size_type nonZeros )
206 const size_type needed = values_.size() + nonZeros ;
207 if( values_.capacity() < needed )
209 const size_type estimate = needed * 1.1;
210 values_.reserve( estimate );
211 cols_.reserve( estimate );
215 void push_back(
const block_type& value,
const size_type index )
217 values_.push_back( value );
218 cols_.push_back( index );
221 std::vector< size_type > rows_;
222 std::vector< block_type > values_;
223 std::vector< size_type > cols_;
231 category = std::is_same<ParallelInfoT, Dune::Amg::SequentialInformation>::value ?
232 Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping
242 template<
class BlockType,
class Alloc>
248 comm_(nullptr),
w_(w),
249 relaxation_( std::abs( w - 1.0 ) > 1e-15 )
253 init( reinterpret_cast<const Matrix&>(A), n );
262 template<
class BlockType,
class Alloc>
276 template<
class BlockType,
class Alloc>
283 relaxation_( std::abs( w - 1.0 ) > 1e-15 )
287 init( reinterpret_cast<const Matrix&>(A), 0 );
295 virtual void pre (Domain& x, Range& b)
297 DUNE_UNUSED_PARAMETER(x);
298 DUNE_UNUSED_PARAMETER(b);
306 virtual void apply (Domain& v,
const Range& d)
308 Range& md =
const_cast<Range&
>(d);
309 copyOwnerToAll( md );
312 typedef typename Range ::block_type dblock;
313 typedef typename Domain::block_type vblock;
315 const size_type iEnd =
lower_.rows();
316 const size_type lastRow = iEnd - 1;
317 if( iEnd != upper_.rows() )
324 for( size_type i=0; i<iEnd; ++ i )
326 dblock rhs( d[ i ] );
327 const size_type rowI =
lower_.rows_[ i ];
328 const size_type rowINext =
lower_.rows_[ i+1 ];
330 for( size_type col = rowI; col < rowINext; ++ col )
332 lower_.values_[ col ].mmv( v[
lower_.cols_[ col ] ], rhs );
340 for( size_type i=0; i<iEnd; ++ i )
342 vblock& vBlock = v[ lastRow - i ];
343 vblock rhs ( vBlock );
344 const size_type rowI = upper_.rows_[ i ];
345 const size_type rowINext = upper_.rows_[ i+1 ];
347 for( size_type col = rowI; col < rowINext; ++ col )
349 upper_.values_[ col ].mmv( v[ upper_.cols_[ col ] ], rhs );
353 inv_[ i ].mv( rhs, vBlock);
364 void copyOwnerToAll( V& v )
const
367 comm_->copyOwnerToAll(v, v);
378 DUNE_UNUSED_PARAMETER(x);
382 void init(
const Matrix& A,
const int iluIteration )
384 int ilu_setup_successful = 1;
386 const int rank = ( comm_ ) ? comm_->communicator().rank() : 0;
388 std::unique_ptr< Matrix > ILU;
392 if( iluIteration == 0 ) {
394 ILU.reset(
new Matrix( A ) );
395 bilu0_decomposition( *ILU );
399 ILU.reset(
new Matrix( A.N(), A.M(), Matrix::row_wise) );
400 bilu_decomposition( A, iluIteration, *ILU );
403 catch ( Dune::MatrixBlockError error )
405 message = error.what();
406 std::cerr<<
"Exception occured on process " << rank <<
" during " <<
407 "setup of ILU0 preconditioner with message: " <<
409 ilu_setup_successful = 0;
413 if ( comm_ && comm_->communicator().min(ilu_setup_successful) == 0 )
415 throw Dune::MatrixBlockError();
419 detail::convertToCRS( *ILU,
lower_, upper_, inv_ );
426 std::vector< block_type > inv_;
428 const ParallelInfo* comm_;
431 const bool relaxation_;
virtual void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: ParallelOverlappingILU0.hpp:295
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:35
Range range_type
The range type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:175
const field_type w_
The relaxation factor to use.
Definition: ParallelOverlappingILU0.hpp:430
CRS lower_
The ILU0 decomposition of the matrix.
Definition: ParallelOverlappingILU0.hpp:424
Domain domain_type
The domain type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:173
std::remove_const< Matrix >::type matrix_type
The matrix type the preconditioner is for.
Definition: ParallelOverlappingILU0.hpp:171
The category the preconditioner is part of.
Definition: ParallelOverlappingILU0.hpp:231
virtual void apply(Domain &v, const Range &d)
Apply the preconditoner.
Definition: ParallelOverlappingILU0.hpp:306
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const int n, const field_type w)
Constructor.
Definition: ParallelOverlappingILU0.hpp:243
Domain::field_type field_type
The field type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:177
Definition: ParallelOverlappingILU0.hpp:183
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const field_type w)
Constructor.
Definition: ParallelOverlappingILU0.hpp:263
virtual void post(Range &x)
Clean up.
Definition: ParallelOverlappingILU0.hpp:376
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const ParallelInfo &comm, const field_type w)
Constructor.
Definition: ParallelOverlappingILU0.hpp:277