20 #ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
21 #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
23 #include <opm/simulators/linalg/GraphColoring.hpp>
24 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
25 #include <opm/common/ErrorMacros.hpp>
26 #include <dune/common/version.hh>
27 #include <dune/istl/preconditioner.hh>
28 #include <dune/istl/ilu.hh>
29 #include <dune/istl/paamg/smoother.hh>
30 #include <dune/istl/paamg/graph.hh>
31 #include <dune/istl/paamg/pinfo.hh>
33 #include <type_traits>
44 template<
class Matrix,
class Domain,
class Range,
class ParallelInfo = Dune::Amg::SequentialInformation>
45 class ParallelOverlappingILU0;
62 if( 0 == milu.compare(
"MILU_1") )
66 if ( 0 == milu.compare(
"MILU_2") )
70 if ( 0 == milu.compare(
"MILU_3") )
79 :
public Dune::Amg::DefaultSmootherArgs<F>
114 template<
class M,
class X,
class Y,
class C>
115 struct SmootherTraits<
Opm::ParallelOverlappingILU0<M,X,Y,C> >
126 template<
class Matrix,
class Domain,
class Range,
class ParallelInfo>
127 struct ConstructionTraits<
Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> >
130 typedef DefaultParallelConstructionArgs<T,ParallelInfo> Arguments;
132 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
141 new T(args.getMatrix(),
143 args.getArgs().getN(),
144 args.getArgs().relaxationFactor,
145 args.getArgs().getMilu()) );
148 #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
150 static inline void deconstruct(
T* bp)
170 virtual std::size_t operator[](std::size_t i)
const = 0;
176 virtual std::size_t operator[](std::size_t i)
const
185 : ordering_(&ordering)
187 virtual std::size_t operator[](std::size_t i)
const
189 return (*ordering_)[i];
191 const std::vector<std::size_t>* ordering_;
197 T operator()(
const T& t)
206 T operator()(
const T&)
214 double operator()(
const T& t)
230 double operator()(
const T& t)
245 T operator()(
const T& t)
252 template<
class M,
class F1=detail::IdentityFunctor,
class F2=detail::OneFunctor >
253 void milu0_decomposition(M& A, F1 absFunctor = F1(), F2 signFunctor = F2(),
254 std::vector<typename M::block_type>* diagonal =
nullptr)
258 diagonal->reserve(A.N());
261 for (
auto irow = A.begin(), iend = A.end(); irow != iend; ++irow)
263 auto a_i_end = irow->end();
264 auto a_ik = irow->begin();
266 std::array<typename M::field_type, M::block_type::rows> sum_dropped{};
270 for ( ; a_ik.index() < irow.index(); ++a_ik )
272 auto k = a_ik.index();
273 auto a_kk = A[k].find(k);
275 a_ik->rightmultiply(*a_kk);
279 auto a_k_end = A[k].end();
280 auto a_kj = a_kk, a_ij = a_ik;
283 while ( a_kj != a_k_end)
285 auto modifier = *a_kj;
286 modifier.leftmultiply(*a_ik);
288 while( a_ij != a_i_end && a_ij.index() < a_kj.index())
293 if ( a_ij != a_i_end && a_ij.index() == a_kj.index() )
301 auto entry = sum_dropped.begin();
302 for(
const auto& row: modifier )
304 for(
const auto& colEntry: row )
306 *entry += absFunctor(-colEntry);
315 if ( a_ik.index() != irow.index() )
316 OPM_THROW(std::logic_error,
"Matrix is missing diagonal for row " << irow.index());
319 for(
const auto& entry: sum_dropped)
321 auto& bdiag = (*a_ik)[index][index];
322 bdiag += signFunctor(bdiag) * entry;
328 diagonal->push_back(*a_ik);
335 void milu0_decomposition(M& A,
336 std::vector<typename M::block_type>* diagonal)
338 milu0_decomposition(A, detail::IdentityFunctor(), detail::OneFunctor(),
343 void milun_decomposition(
const M& A,
int n,
MILU_VARIANT milu, M& ILU,
344 Reorderer& ordering, Reorderer& inverseOrdering)
346 using Map = std::map<std::size_t, int>;
348 auto iluRow =
ILU.createbegin();
350 for(std::size_t i = 0, iend = A.N(); i < iend; ++i)
352 auto& orow = A[inverseOrdering[i]];
355 for (
auto col = orow.begin(), cend = orow.end(); col != cend; ++col)
357 rowPattern[ordering[col.index()]] = 0;
360 for(
auto ik = rowPattern.begin(); ik->first < i; ++ik)
362 if ( ik->second < n )
364 auto& rowk =
ILU[ik->first];
366 for (
auto kj = rowk.find(ik->first), endk = rowk.end();
371 int generation = (*kj)[0][0];
374 auto ij = rowPattern.find(kj.index());
375 if ( ij == rowPattern.end() )
377 rowPattern[ordering[kj.index()]] = generation + 1;
384 for (
const auto& entry : rowPattern)
386 iluRow.insert(entry.first);
391 auto generationPair = rowPattern.begin();
392 for (
auto col = ILU[i].begin(), cend = ILU[i].end(); col != cend;
393 ++col, ++generationPair)
395 assert(col.index() == generationPair->first);
396 (*col)[0][0] = generationPair->second;
401 for(
auto iter=A.begin(), iend = A.end(); iter != iend; ++iter)
403 auto& newRow =
ILU[ordering[iter.index()]];
405 for (
auto& col: newRow)
410 for(
auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
412 newRow[ordering[col.index()]] = *col;
419 detail::milu0_decomposition ( ILU);
422 detail::milu0_decomposition ( ILU, detail::IdentityFunctor(),
423 detail::SignFunctor() );
426 detail::milu0_decomposition ( ILU, detail::AbsFunctor(),
427 detail::SignFunctor() );
430 detail::milu0_decomposition ( ILU, detail::IdentityFunctor(),
431 detail::IsPositiveFunctor() );
434 #if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
435 bilu0_decomposition( ILU );
437 Dune::ILU::blockILU0Decomposition( ILU );
445 void ghost_last_bilu0_decomposition (M& A,
size_t interiorSize)
448 typedef typename M::RowIterator rowiterator;
449 typedef typename M::ColIterator coliterator;
450 typedef typename M::block_type block;
453 for (rowiterator i = A.begin(); i.index() < interiorSize; ++i)
456 coliterator endij=(*i).end();
460 for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
463 coliterator jj = A[ij.index()].find(ij.index());
466 (*ij).rightmultiply(*jj);
469 coliterator endjk=A[ij.index()].end();
470 coliterator jk=jj; ++jk;
471 coliterator ik=ij; ++ik;
472 while (ik!=endij && jk!=endjk)
473 if (ik.index()==jk.index())
482 if (ik.index()<jk.index())
490 if (ij.index()!=i.index())
491 DUNE_THROW(Dune::ISTLError,
"diagonal entry missing");
495 catch (Dune::FMatrixError & e) {
496 DUNE_THROW(Dune::ISTLError,
"ILU failed to invert matrix block");
502 template<
class M,
class CRS,
class InvVector>
503 void convertToCRS(
const M& A, CRS& lower, CRS& upper, InvVector& inv )
512 typedef typename M :: size_type size_type;
517 lower.resize( A.N() );
518 upper.resize( A.N() );
522 size_type numLower = 0;
523 size_type numUpper = 0;
524 const auto endi = A.end();
525 for (
auto i = A.begin(); i != endi; ++i) {
526 const size_type iIndex = i.index();
527 size_type numLowerRow = 0;
528 for (
auto j = (*i).begin(); j.index() < iIndex; ++j) {
531 numLower += numLowerRow;
532 numUpper += (*i).size() - numLowerRow - 1;
534 assert(numLower + numUpper + A.N() == A.nonzeroes());
536 lower.reserveAdditional( numLower );
540 size_type colcount = 0;
541 lower.rows_[ 0 ] = colcount;
542 for (
auto i=A.begin(); i!=endi; ++i, ++row)
544 const size_type iIndex = i.index();
547 for (
auto j=(*i).begin(); j.index() < iIndex; ++j )
549 lower.push_back( (*j), j.index() );
552 lower.rows_[ iIndex+1 ] = colcount;
555 assert(colcount == numLower);
557 const auto rendi = A.beforeBegin();
560 upper.rows_[ 0 ] = colcount ;
562 upper.reserveAdditional( numUpper );
566 for (
auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
568 const size_type iIndex = i.index();
572 for (
auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
574 const size_type jIndex = j.index();
575 if( j.index() == iIndex )
580 else if ( j.index() >= i.index() )
582 upper.push_back( (*j), jIndex );
586 upper.rows_[ row+1 ] = colcount;
588 assert(colcount == numUpper);
608 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
612 typedef ParallelInfoT ParallelInfo;
625 typedef typename matrix_type::block_type block_type;
626 typedef typename matrix_type::size_type size_type;
631 CRS() : nRows_( 0 ) {}
633 size_type rows()
const {
return nRows_; }
635 size_type nonZeros()
const
637 assert( rows_[ rows() ] != size_type(-1) );
638 return rows_[ rows() ];
641 void resize(
const size_type nRows )
643 if( nRows_ != nRows )
646 rows_.resize( nRows_+1, size_type(-1) );
650 void reserveAdditional(
const size_type nonZeros )
652 const size_type needed = values_.size() + nonZeros ;
653 if( values_.capacity() < needed )
655 const size_type estimate = needed * 1.1;
656 values_.reserve( estimate );
657 cols_.reserve( estimate );
661 void push_back(
const block_type& value,
const size_type index )
663 values_.push_back( value );
664 cols_.push_back( index );
675 std::vector< size_type > rows_;
676 std::vector< block_type > values_;
677 std::vector< size_type > cols_;
682 Dune::SolverCategory::Category category()
const override
684 return std::is_same<ParallelInfoT, Dune::Amg::SequentialInformation>::value ?
685 Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
701 template<
class BlockType,
class Alloc>
705 bool reorder_sphere=
true)
709 comm_(nullptr),
w_(w),
710 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
711 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
712 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
714 interiorSize_ = A.N();
732 template<
class BlockType,
class Alloc>
734 const ParallelInfo& comm,
const int n,
const field_type w,
736 bool reorder_sphere=
true)
741 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
742 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
743 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
745 interiorSize_ = A.N();
763 template<
class BlockType,
class Alloc>
766 bool reorder_sphere=
true)
784 template<
class BlockType,
class Alloc>
788 bool reorder_sphere=
true)
793 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
794 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
795 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
797 interiorSize_ = A.N();
817 template<
class BlockType,
class Alloc>
819 const ParallelInfo& comm,
821 size_type interiorSize,
bool redblack=
false,
822 bool reorder_sphere=
true)
827 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
828 interiorSize_(interiorSize),
829 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
830 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
842 virtual void pre (Domain& x, Range& b)
override
844 DUNE_UNUSED_PARAMETER(x);
845 DUNE_UNUSED_PARAMETER(b);
853 virtual void apply (Domain& v,
const Range& d)
override
859 typedef typename Range ::block_type dblock;
860 typedef typename Domain::block_type vblock;
862 const size_type iEnd =
lower_.rows();
863 const size_type lastRow = iEnd - 1;
864 size_type upperLoppStart = iEnd - interiorSize_;
865 size_type lowerLoopEnd = interiorSize_;
866 if( iEnd != upper_.rows() )
868 OPM_THROW(std::logic_error,
"ILU: number of lower and upper rows must be the same");
872 for( size_type i=0; i<lowerLoopEnd; ++ i )
874 dblock rhs( md[ i ] );
875 const size_type rowI =
lower_.rows_[ i ];
876 const size_type rowINext =
lower_.rows_[ i+1 ];
878 for( size_type col = rowI; col < rowINext; ++ col )
880 lower_.values_[ col ].mmv( mv[
lower_.cols_[ col ] ], rhs );
886 for( size_type i=upperLoppStart; i<iEnd; ++ i )
888 vblock& vBlock = mv[ lastRow - i ];
889 vblock rhs ( vBlock );
890 const size_type rowI = upper_.rows_[ i ];
891 const size_type rowINext = upper_.rows_[ i+1 ];
893 for( size_type col = rowI; col < rowINext; ++ col )
895 upper_.values_[ col ].mmv( mv[ upper_.cols_[ col ] ], rhs );
899 inv_[ i ].mv( rhs, vBlock);
902 copyOwnerToAll( mv );
911 void copyOwnerToAll( V& v )
const
914 comm_->copyOwnerToAll(v, v);
923 virtual void post (Range& x)
override
925 DUNE_UNUSED_PARAMETER(x);
928 virtual void update()
override
933 if ( comm_ && comm_->communicator().size()<=0 )
937 OPM_THROW(std::logic_error,
"Expected a matrix with zero rows for an invalid communicator.");
946 int ilu_setup_successful = 1;
948 const int rank = ( comm_ ) ? comm_->communicator().rank() : 0;
950 std::unique_ptr< Matrix >
ILU;
954 using Graph = Dune::Amg::MatrixGraph<const Matrix>;
957 const auto& colors = std::get<0>(colorsTuple);
958 const auto& verticesPerColor = std::get<2>(colorsTuple);
959 auto noColors = std::get<1>(colorsTuple);
960 if ( reorderSphere_ )
972 std::vector<std::size_t> inverseOrdering(
ordering_.size());
973 std::size_t index = 0;
976 inverseOrdering[newIndex] = index++;
981 if( iluIteration_ == 0 ) {
985 ILU.reset(
new Matrix( *A_ ) );
989 ILU.reset(
new Matrix(A_->N(), A_->M(), A_->nonzeroes(), Matrix::row_wise));
992 for(
auto iter=newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
994 const auto& row = (*A_)[inverseOrdering[iter.index()]];
995 for(
auto col = row.begin(), cend = row.end(); col != cend; ++col)
1001 for(
auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
1003 auto& newRow = newA[
ordering_[iter.index()]];
1004 for(
auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
1014 detail::milu0_decomposition ( *ILU);
1017 detail::milu0_decomposition ( *ILU, detail::IdentityFunctor(),
1018 detail::SignFunctor() );
1021 detail::milu0_decomposition ( *ILU, detail::AbsFunctor(),
1022 detail::SignFunctor() );
1025 detail::milu0_decomposition ( *ILU, detail::IdentityFunctor(),
1026 detail::IsPositiveFunctor() );
1029 if (interiorSize_ == A_->N())
1030 #if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
1031 bilu0_decomposition( *ILU );
1033 Dune::ILU::blockILU0Decomposition( *ILU );
1036 detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
1042 ILU.reset(
new Matrix( A_->N(), A_->M(), Matrix::row_wise) );
1043 std::unique_ptr<detail::Reorderer> reorderer, inverseReorderer;
1046 reorderer.reset(
new detail::NoReorderer());
1047 inverseReorderer.reset(
new detail::NoReorderer());
1051 reorderer.reset(
new detail::RealReorderer(
ordering_));
1052 inverseReorderer.reset(
new detail::RealReorderer(inverseOrdering));
1055 milun_decomposition( *A_, iluIteration_, milu_, *ILU, *reorderer, *inverseReorderer );
1058 catch (
const Dune::MatrixBlockError& error)
1060 message = error.what();
1061 std::cerr<<
"Exception occurred on process " << rank <<
" during " <<
1062 "setup of ILU0 preconditioner with message: " <<
1064 ilu_setup_successful = 0;
1068 const bool parallel_failure = comm_ && comm_->communicator().min(ilu_setup_successful) == 0;
1069 const bool local_failure = ilu_setup_successful == 0;
1070 if ( local_failure || parallel_failure )
1072 throw Dune::MatrixBlockError();
1076 detail::convertToCRS( *ILU,
lower_, upper_, inv_ );
1089 return const_cast<Range&
>(d);
1122 void reorderBack(
const Range& reorderedV, Range& v)
1129 v[i++] = reorderedV[index];
1137 std::vector< block_type > inv_;
1145 const ParallelInfo* comm_;
1148 const bool relaxation_;
1149 size_type interiorSize_;
1154 bool reorderSphere_;
Interface class adding the update() method to the preconditioner interface.
Definition: PreconditionerWithUpdate.hpp:32
Definition: ParallelOverlappingILU0.hpp:80
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:611
virtual void post(Range &x) override
Clean up.
Definition: ParallelOverlappingILU0.hpp:923
Range & reorderD(const Range &d)
Reorder D if needed and return a reference to it.
Definition: ParallelOverlappingILU0.hpp:1081
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const field_type w, MILU_VARIANT milu, bool redblack=false, bool reorder_sphere=true)
Constructor.
Definition: ParallelOverlappingILU0.hpp:764
Domain & reorderV(Domain &v)
Reorder V if needed and return a reference to it.
Definition: ParallelOverlappingILU0.hpp:1104
std::vector< std::size_t > ordering_
the reordering of the unknowns
Definition: ParallelOverlappingILU0.hpp:1139
std::remove_const< Matrix >::type matrix_type
The matrix type the preconditioner is for.
Definition: ParallelOverlappingILU0.hpp:617
Domain reorderedV_
The reordered left hand side.
Definition: ParallelOverlappingILU0.hpp:1143
Range reorderedD_
The reordered right hand side.
Definition: ParallelOverlappingILU0.hpp:1141
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const ParallelInfo &comm, const field_type w, MILU_VARIANT milu, bool redblack=false, bool reorder_sphere=true)
Constructor.
Definition: ParallelOverlappingILU0.hpp:785
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const ParallelInfo &comm, const field_type w, MILU_VARIANT milu, size_type interiorSize, bool redblack=false, bool reorder_sphere=true)
Constructor.
Definition: ParallelOverlappingILU0.hpp:818
const field_type w_
The relaxation factor to use.
Definition: ParallelOverlappingILU0.hpp:1147
virtual void apply(Domain &v, const Range &d) override
Apply the preconditoner.
Definition: ParallelOverlappingILU0.hpp:853
Domain::field_type field_type
The field type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:623
CRS lower_
The ILU0 decomposition of the matrix.
Definition: ParallelOverlappingILU0.hpp:1135
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const int n, const field_type w, MILU_VARIANT milu, bool redblack=false, bool reorder_sphere=true)
Constructor.
Definition: ParallelOverlappingILU0.hpp:702
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const ParallelInfo &comm, const int n, const field_type w, MILU_VARIANT milu, bool redblack=false, bool reorder_sphere=true)
Constructor gets all parameters to operate the prec.
Definition: ParallelOverlappingILU0.hpp:733
Domain domain_type
The domain type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:619
virtual void pre(Domain &x, Range &b) override
Prepare the preconditioner.
Definition: ParallelOverlappingILU0.hpp:842
Range range_type
The range type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:621
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:26
MILU_VARIANT
Definition: ParallelOverlappingILU0.hpp:47
@ MILU_1
sum(dropped entries)
@ MILU_2
sum(dropped entries)
@ MILU_3
sum(|dropped entries|)
@ MILU_4
sum(dropped entries)
@ ILU
Do not perform modified ILU.
std::vector< std::size_t > reorderVerticesPreserving(const std::vector< int > &colors, int noColors, const std::vector< std::size_t > &verticesPerColor, const Graph &graph)
! Reorder colored graph preserving order of vertices with the same color.
Definition: GraphColoring.hpp:186
std::vector< std::size_t > reorderVerticesSpheres(const std::vector< int > &colors, int noColors, const std::vector< std::size_t > &verticesPerColor, const Graph &graph, typename Graph::VertexDescriptor root)
! Reorder Vetrices in spheres
Definition: GraphColoring.hpp:206
std::tuple< std::vector< int >, int, std::vector< std::size_t > > colorVerticesWelshPowell(const Graph &graph)
Color the vertices of graph.
Definition: GraphColoring.hpp:118
Definition: ParallelOverlappingILU0.hpp:630
Definition: ParallelOverlappingILU0.hpp:243
Definition: ParallelOverlappingILU0.hpp:195
Definition: ParallelOverlappingILU0.hpp:228
Definition: ParallelOverlappingILU0.hpp:175
Definition: ParallelOverlappingILU0.hpp:204
Definition: ParallelOverlappingILU0.hpp:183
Definition: ParallelOverlappingILU0.hpp:169
Definition: ParallelOverlappingILU0.hpp:212