My Project
ParallelOverlappingILU0.hpp
1 /*
2  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
3  Copyright 2015 Statoil AS
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 #ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
21 #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
22 
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>
32 
33 #include <type_traits>
34 #include <numeric>
35 #include <limits>
36 #include <cstddef>
37 #include <string>
38 
39 namespace Opm
40 {
41 
42 //template<class M, class X, class Y, class C>
43 //class ParallelOverlappingILU0;
44 template<class Matrix, class Domain, class Range, class ParallelInfo = Dune::Amg::SequentialInformation>
45 class ParallelOverlappingILU0;
46 
47 enum class MILU_VARIANT{
49  ILU = 0,
51  MILU_1 = 1,
53  MILU_2 = 2,
55  MILU_3 = 3,
57  MILU_4 = 4
58 };
59 
60 inline MILU_VARIANT convertString2Milu(std::string milu)
61 {
62  if( 0 == milu.compare("MILU_1") )
63  {
64  return MILU_VARIANT::MILU_1;
65  }
66  if ( 0 == milu.compare("MILU_2") )
67  {
68  return MILU_VARIANT::MILU_2;
69  }
70  if ( 0 == milu.compare("MILU_3") )
71  {
72  return MILU_VARIANT::MILU_3;
73  }
74  return MILU_VARIANT::ILU;
75 }
76 
77 template<class F>
79  : public Dune::Amg::DefaultSmootherArgs<F>
80 {
81  public:
83  : milu_(milu), n_(0)
84  {}
85  void setMilu(MILU_VARIANT milu)
86  {
87  milu_ = milu;
88  }
89  MILU_VARIANT getMilu() const
90  {
91  return milu_;
92  }
93  void setN(int n)
94  {
95  n_ = n;
96  }
97  int getN() const
98  {
99  return n_;
100  }
101  private:
102  MILU_VARIANT milu_;
103  int n_;
104 };
105 } // end namespace Opm
106 
107 namespace Dune
108 {
109 
110 namespace Amg
111 {
112 
113 
114 template<class M, class X, class Y, class C>
115 struct SmootherTraits<Opm::ParallelOverlappingILU0<M,X,Y,C> >
116 {
118 };
119 
126 template<class Matrix, class Domain, class Range, class ParallelInfo>
127 struct ConstructionTraits<Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> >
128 {
130  typedef DefaultParallelConstructionArgs<T,ParallelInfo> Arguments;
131 
132 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
133  typedef std::shared_ptr< T > ParallelOverlappingILU0Pointer;
134 #else
136 #endif
137 
138  static inline ParallelOverlappingILU0Pointer construct(Arguments& args)
139  {
141  new T(args.getMatrix(),
142  args.getComm(),
143  args.getArgs().getN(),
144  args.getArgs().relaxationFactor,
145  args.getArgs().getMilu()) );
146  }
147 
148 #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
149  // this method is not needed anymore in 2.7 since std::shared_ptr is used
150  static inline void deconstruct(T* bp)
151  {
152  delete bp;
153  }
154 #endif
155 
156 };
157 
158 } // end namespace Amg
159 
160 
161 
162 } // end namespace Dune
163 
164 namespace Opm
165 {
166  namespace detail
167  {
168  struct Reorderer
169  {
170  virtual std::size_t operator[](std::size_t i) const = 0;
171  virtual ~Reorderer() {}
172  };
173 
174  struct NoReorderer : public Reorderer
175  {
176  virtual std::size_t operator[](std::size_t i) const
177  {
178  return i;
179  }
180  };
181 
182  struct RealReorderer : public Reorderer
183  {
184  RealReorderer(const std::vector<std::size_t>& ordering)
185  : ordering_(&ordering)
186  {}
187  virtual std::size_t operator[](std::size_t i) const
188  {
189  return (*ordering_)[i];
190  }
191  const std::vector<std::size_t>* ordering_;
192  };
193 
195  {
196  template<class T>
197  T operator()(const T& t)
198  {
199  return t;
200  }
201  };
202 
203  struct OneFunctor
204  {
205  template<class T>
206  T operator()(const T&)
207  {
208  return 1.0;
209  }
210  };
211  struct SignFunctor
212  {
213  template<class T>
214  double operator()(const T& t)
215  {
216  if ( t < 0.0 )
217  {
218  return -1;
219  }
220  else
221  {
222  return 1.0;
223  }
224  }
225  };
226 
228  {
229  template<class T>
230  double operator()(const T& t)
231  {
232  if ( t < 0.0 )
233  {
234  return 0;
235  }
236  else
237  {
238  return 1;
239  }
240  }
241  };
242  struct AbsFunctor
243  {
244  template<class T>
245  T operator()(const T& t)
246  {
247  using std::abs;
248  return abs(t);
249  }
250  };
251 
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)
255  {
256  if( diagonal )
257  {
258  diagonal->reserve(A.N());
259  }
260 
261  for ( auto irow = A.begin(), iend = A.end(); irow != iend; ++irow)
262  {
263  auto a_i_end = irow->end();
264  auto a_ik = irow->begin();
265 
266  std::array<typename M::field_type, M::block_type::rows> sum_dropped{};
267 
268  // Eliminate entries in lower triangular matrix
269  // and store factors for L
270  for ( ; a_ik.index() < irow.index(); ++a_ik )
271  {
272  auto k = a_ik.index();
273  auto a_kk = A[k].find(k);
274  // L_ik = A_kk^-1 * A_ik
275  a_ik->rightmultiply(*a_kk);
276 
277  // modify the rest of the row, everything right of a_ik
278  // a_i* -=a_ik * a_k*
279  auto a_k_end = A[k].end();
280  auto a_kj = a_kk, a_ij = a_ik;
281  ++a_kj; ++a_ij;
282 
283  while ( a_kj != a_k_end)
284  {
285  auto modifier = *a_kj;
286  modifier.leftmultiply(*a_ik);
287 
288  while( a_ij != a_i_end && a_ij.index() < a_kj.index())
289  {
290  ++a_ij;
291  }
292 
293  if ( a_ij != a_i_end && a_ij.index() == a_kj.index() )
294  {
295  // Value is not dropped
296  *a_ij -= modifier;
297  ++a_ij; ++a_kj;
298  }
299  else
300  {
301  auto entry = sum_dropped.begin();
302  for( const auto& row: modifier )
303  {
304  for( const auto& colEntry: row )
305  {
306  *entry += absFunctor(-colEntry);
307  }
308  ++entry;
309  }
310  ++a_kj;
311  }
312  }
313  }
314 
315  if ( a_ik.index() != irow.index() )
316  OPM_THROW(std::logic_error, "Matrix is missing diagonal for row " << irow.index());
317 
318  int index = 0;
319  for(const auto& entry: sum_dropped)
320  {
321  auto& bdiag = (*a_ik)[index][index];
322  bdiag += signFunctor(bdiag) * entry;
323  ++index;
324  }
325 
326  if ( diagonal )
327  {
328  diagonal->push_back(*a_ik);
329  }
330  a_ik->invert(); // compute inverse of diagonal block
331  }
332  }
333 
334  template<class M>
335  void milu0_decomposition(M& A,
336  std::vector<typename M::block_type>* diagonal)
337  {
338  milu0_decomposition(A, detail::IdentityFunctor(), detail::OneFunctor(),
339  diagonal);
340  }
341 
342  template<class M>
343  void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU,
344  Reorderer& ordering, Reorderer& inverseOrdering)
345  {
346  using Map = std::map<std::size_t, int>;
347 
348  auto iluRow = ILU.createbegin();
349 
350  for(std::size_t i = 0, iend = A.N(); i < iend; ++i)
351  {
352  auto& orow = A[inverseOrdering[i]];
353 
354  Map rowPattern;
355  for ( auto col = orow.begin(), cend = orow.end(); col != cend; ++col)
356  {
357  rowPattern[ordering[col.index()]] = 0;
358  }
359 
360  for(auto ik = rowPattern.begin(); ik->first < i; ++ik)
361  {
362  if ( ik->second < n )
363  {
364  auto& rowk = ILU[ik->first];
365 
366  for ( auto kj = rowk.find(ik->first), endk = rowk.end();
367  kj != endk; ++kj)
368  {
369  // Assume double and block_type FieldMatrix
370  // first element is misused to store generation number
371  int generation = (*kj)[0][0];
372  if(generation < n)
373  {
374  auto ij = rowPattern.find(kj.index());
375  if ( ij == rowPattern.end() )
376  {
377  rowPattern[ordering[kj.index()]] = generation + 1;
378  }
379  }
380  }
381  }
382  }
383  // create the row
384  for (const auto& entry : rowPattern)
385  {
386  iluRow.insert(entry.first);
387  }
388  ++iluRow;
389 
390  // write generation to newly created row.
391  auto generationPair = rowPattern.begin();
392  for ( auto col = ILU[i].begin(), cend = ILU[i].end(); col != cend;
393  ++col, ++generationPair)
394  {
395  assert(col.index() == generationPair->first);
396  (*col)[0][0] = generationPair->second;
397  }
398  }
399 
400  // copy Entries from A
401  for(auto iter=A.begin(), iend = A.end(); iter != iend; ++iter)
402  {
403  auto& newRow = ILU[ordering[iter.index()]];
404  // reset stored generation
405  for ( auto& col: newRow)
406  {
407  col = 0;
408  }
409  // copy row.
410  for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
411  {
412  newRow[ordering[col.index()]] = *col;
413  }
414  }
415  // call decomposition on pattern
416  switch ( milu )
417  {
419  detail::milu0_decomposition ( ILU);
420  break;
422  detail::milu0_decomposition ( ILU, detail::IdentityFunctor(),
423  detail::SignFunctor() );
424  break;
426  detail::milu0_decomposition ( ILU, detail::AbsFunctor(),
427  detail::SignFunctor() );
428  break;
430  detail::milu0_decomposition ( ILU, detail::IdentityFunctor(),
431  detail::IsPositiveFunctor() );
432  break;
433  default:
434 #if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
435  bilu0_decomposition( ILU );
436 #else
437  Dune::ILU::blockILU0Decomposition( ILU );
438 #endif
439  break;
440  }
441  }
442 
444  template<class M>
445  void ghost_last_bilu0_decomposition (M& A, size_t interiorSize)
446  {
447  // iterator types
448  typedef typename M::RowIterator rowiterator;
449  typedef typename M::ColIterator coliterator;
450  typedef typename M::block_type block;
451 
452  // implement left looking variant with stored inverse
453  for (rowiterator i = A.begin(); i.index() < interiorSize; ++i)
454  {
455  // coliterator is diagonal after the following loop
456  coliterator endij=(*i).end(); // end of row i
457  coliterator ij;
458 
459  // eliminate entries left of diagonal; store L factor
460  for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
461  {
462  // find A_jj which eliminates A_ij
463  coliterator jj = A[ij.index()].find(ij.index());
464 
465  // compute L_ij = A_jj^-1 * A_ij
466  (*ij).rightmultiply(*jj);
467 
468  // modify row
469  coliterator endjk=A[ij.index()].end(); // end of row j
470  coliterator jk=jj; ++jk;
471  coliterator ik=ij; ++ik;
472  while (ik!=endij && jk!=endjk)
473  if (ik.index()==jk.index())
474  {
475  block B(*jk);
476  B.leftmultiply(*ij);
477  *ik -= B;
478  ++ik; ++jk;
479  }
480  else
481  {
482  if (ik.index()<jk.index())
483  ++ik;
484  else
485  ++jk;
486  }
487  }
488 
489  // invert pivot and store it in A
490  if (ij.index()!=i.index())
491  DUNE_THROW(Dune::ISTLError,"diagonal entry missing");
492  try {
493  (*ij).invert(); // compute inverse of diagonal block
494  }
495  catch (Dune::FMatrixError & e) {
496  DUNE_THROW(Dune::ISTLError,"ILU failed to invert matrix block");
497  }
498  }
499  }
500 
502  template<class M, class CRS, class InvVector>
503  void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv )
504  {
505  // No need to do anything for 0 rows. Return to prevent indexing a
506  // a zero sized array.
507  if ( A.N() == 0 )
508  {
509  return;
510  }
511 
512  typedef typename M :: size_type size_type;
513 
514  lower.clear();
515  upper.clear();
516  inv.clear();
517  lower.resize( A.N() );
518  upper.resize( A.N() );
519  inv.resize( A.N() );
520 
521  // Count the lower and upper matrix entries.
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) {
529  ++numLowerRow;
530  }
531  numLower += numLowerRow;
532  numUpper += (*i).size() - numLowerRow - 1;
533  }
534  assert(numLower + numUpper + A.N() == A.nonzeroes());
535 
536  lower.reserveAdditional( numLower );
537 
538  // implement left looking variant with stored inverse
539  size_type row = 0;
540  size_type colcount = 0;
541  lower.rows_[ 0 ] = colcount;
542  for (auto i=A.begin(); i!=endi; ++i, ++row)
543  {
544  const size_type iIndex = i.index();
545 
546  // eliminate entries left of diagonal; store L factor
547  for (auto j=(*i).begin(); j.index() < iIndex; ++j )
548  {
549  lower.push_back( (*j), j.index() );
550  ++colcount;
551  }
552  lower.rows_[ iIndex+1 ] = colcount;
553  }
554 
555  assert(colcount == numLower);
556 
557  const auto rendi = A.beforeBegin();
558  row = 0;
559  colcount = 0;
560  upper.rows_[ 0 ] = colcount ;
561 
562  upper.reserveAdditional( numUpper );
563 
564  // NOTE: upper and inv store entries in reverse order, reverse here
565  // relative to ILU
566  for (auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
567  {
568  const size_type iIndex = i.index();
569 
570  // store in reverse row order
571  // eliminate entries left of diagonal; store L factor
572  for (auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
573  {
574  const size_type jIndex = j.index();
575  if( j.index() == iIndex )
576  {
577  inv[ row ] = (*j);
578  break;
579  }
580  else if ( j.index() >= i.index() )
581  {
582  upper.push_back( (*j), jIndex );
583  ++colcount ;
584  }
585  }
586  upper.rows_[ row+1 ] = colcount;
587  }
588  assert(colcount == numUpper);
589  }
590  } // end namespace detail
591 
592 
608 template<class Matrix, class Domain, class Range, class ParallelInfoT>
610  : public Dune::PreconditionerWithUpdate<Domain,Range>
611 {
612  typedef ParallelInfoT ParallelInfo;
613 
614 
615 public:
617  typedef typename std::remove_const<Matrix>::type matrix_type;
619  typedef Domain domain_type;
621  typedef Range range_type;
623  typedef typename Domain::field_type field_type;
624 
625  typedef typename matrix_type::block_type block_type;
626  typedef typename matrix_type::size_type size_type;
627 
628 protected:
629  struct CRS
630  {
631  CRS() : nRows_( 0 ) {}
632 
633  size_type rows() const { return nRows_; }
634 
635  size_type nonZeros() const
636  {
637  assert( rows_[ rows() ] != size_type(-1) );
638  return rows_[ rows() ];
639  }
640 
641  void resize( const size_type nRows )
642  {
643  if( nRows_ != nRows )
644  {
645  nRows_ = nRows ;
646  rows_.resize( nRows_+1, size_type(-1) );
647  }
648  }
649 
650  void reserveAdditional( const size_type nonZeros )
651  {
652  const size_type needed = values_.size() + nonZeros ;
653  if( values_.capacity() < needed )
654  {
655  const size_type estimate = needed * 1.1;
656  values_.reserve( estimate );
657  cols_.reserve( estimate );
658  }
659  }
660 
661  void push_back( const block_type& value, const size_type index )
662  {
663  values_.push_back( value );
664  cols_.push_back( index );
665  }
666 
667  void clear()
668  {
669  rows_.clear();
670  values_.clear();
671  cols_.clear();
672  nRows_= 0;
673  }
674 
675  std::vector< size_type > rows_;
676  std::vector< block_type > values_;
677  std::vector< size_type > cols_;
678  size_type nRows_;
679  };
680 
681 public:
682  Dune::SolverCategory::Category category() const override
683  {
684  return std::is_same<ParallelInfoT, Dune::Amg::SequentialInformation>::value ?
685  Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
686  }
687 
701  template<class BlockType, class Alloc>
702  ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
703  const int n, const field_type w,
704  MILU_VARIANT milu, bool redblack=false,
705  bool reorder_sphere=true)
706  : lower_(),
707  upper_(),
708  inv_(),
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)
713  {
714  interiorSize_ = A.N();
715  // BlockMatrix is a Subclass of FieldMatrix that just adds
716  // methods. Therefore this cast should be safe.
717  update();
718  }
719 
732  template<class BlockType, class Alloc>
733  ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
734  const ParallelInfo& comm, const int n, const field_type w,
735  MILU_VARIANT milu, bool redblack=false,
736  bool reorder_sphere=true)
737  : lower_(),
738  upper_(),
739  inv_(),
740  comm_(&comm), w_(w),
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)
744  {
745  interiorSize_ = A.N();
746  // BlockMatrix is a Subclass of FieldMatrix that just adds
747  // methods. Therefore this cast should be safe.
748  update();
749  }
750 
763  template<class BlockType, class Alloc>
764  ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
765  const field_type w, MILU_VARIANT milu, bool redblack=false,
766  bool reorder_sphere=true)
767  : ParallelOverlappingILU0( A, 0, w, milu, redblack, reorder_sphere )
768  {
769  }
770 
784  template<class BlockType, class Alloc>
785  ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
786  const ParallelInfo& comm, const field_type w,
787  MILU_VARIANT milu, bool redblack=false,
788  bool reorder_sphere=true)
789  : lower_(),
790  upper_(),
791  inv_(),
792  comm_(&comm), w_(w),
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)
796  {
797  interiorSize_ = A.N();
798  // BlockMatrix is a Subclass of FieldMatrix that just adds
799  // methods. Therefore this cast should be safe.
800  update();
801  }
802 
817  template<class BlockType, class Alloc>
818  ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
819  const ParallelInfo& comm,
820  const field_type w, MILU_VARIANT milu,
821  size_type interiorSize, bool redblack=false,
822  bool reorder_sphere=true)
823  : lower_(),
824  upper_(),
825  inv_(),
826  comm_(&comm), w_(w),
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)
831  {
832  // BlockMatrix is a Subclass of FieldMatrix that just adds
833  // methods. Therefore this cast should be safe.
834  update( );
835  }
836 
842  virtual void pre (Domain& x, Range& b) override
843  {
844  DUNE_UNUSED_PARAMETER(x);
845  DUNE_UNUSED_PARAMETER(b);
846  }
847 
853  virtual void apply (Domain& v, const Range& d) override
854  {
855  Range& md = reorderD(d);
856  Domain& mv = reorderV(v);
857 
858  // iterator types
859  typedef typename Range ::block_type dblock;
860  typedef typename Domain::block_type vblock;
861 
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() )
867  {
868  OPM_THROW(std::logic_error,"ILU: number of lower and upper rows must be the same");
869  }
870 
871  // lower triangular solve
872  for( size_type i=0; i<lowerLoopEnd; ++ i )
873  {
874  dblock rhs( md[ i ] );
875  const size_type rowI = lower_.rows_[ i ];
876  const size_type rowINext = lower_.rows_[ i+1 ];
877 
878  for( size_type col = rowI; col < rowINext; ++ col )
879  {
880  lower_.values_[ col ].mmv( mv[ lower_.cols_[ col ] ], rhs );
881  }
882 
883  mv[ i ] = rhs; // Lii = I
884  }
885 
886  for( size_type i=upperLoppStart; i<iEnd; ++ i )
887  {
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 ];
892 
893  for( size_type col = rowI; col < rowINext; ++ col )
894  {
895  upper_.values_[ col ].mmv( mv[ upper_.cols_[ col ] ], rhs );
896  }
897 
898  // apply inverse and store result
899  inv_[ i ].mv( rhs, vBlock);
900  }
901 
902  copyOwnerToAll( mv );
903 
904  if( relaxation_ ) {
905  mv *= w_;
906  }
907  reorderBack(mv, v);
908  }
909 
910  template <class V>
911  void copyOwnerToAll( V& v ) const
912  {
913  if( comm_ ) {
914  comm_->copyOwnerToAll(v, v);
915  }
916  }
917 
923  virtual void post (Range& x) override
924  {
925  DUNE_UNUSED_PARAMETER(x);
926  }
927 
928  virtual void update() override
929  {
930  // (For older DUNE versions the communicator might be
931  // invalid if redistribution in AMG happened on the coarset level.
932  // Therefore we check for nonzero size
933  if ( comm_ && comm_->communicator().size()<=0 )
934  {
935  if ( A_->N() > 0 )
936  {
937  OPM_THROW(std::logic_error, "Expected a matrix with zero rows for an invalid communicator.");
938  }
939  else
940  {
941  // simply set the communicator to null
942  comm_ = nullptr;
943  }
944  }
945 
946  int ilu_setup_successful = 1;
947  std::string message;
948  const int rank = ( comm_ ) ? comm_->communicator().rank() : 0;
949 
950  std::unique_ptr< Matrix > ILU;
951 
952  if ( redBlack_ )
953  {
954  using Graph = Dune::Amg::MatrixGraph<const Matrix>;
955  Graph graph(*A_);
956  auto colorsTuple = colorVerticesWelshPowell(graph);
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_ )
961  {
962  ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor,
963  graph, 0);
964  }
965  else
966  {
967  ordering_ = reorderVerticesPreserving(colors, noColors, verticesPerColor,
968  graph);
969  }
970  }
971 
972  std::vector<std::size_t> inverseOrdering(ordering_.size());
973  std::size_t index = 0;
974  for( auto newIndex: ordering_)
975  {
976  inverseOrdering[newIndex] = index++;
977  }
978 
979  try
980  {
981  if( iluIteration_ == 0 ) {
982  // create ILU-0 decomposition
983  if ( ordering_.empty() )
984  {
985  ILU.reset( new Matrix( *A_ ) );
986  }
987  else
988  {
989  ILU.reset( new Matrix(A_->N(), A_->M(), A_->nonzeroes(), Matrix::row_wise));
990  auto& newA = *ILU;
991  // Create sparsity pattern
992  for(auto iter=newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
993  {
994  const auto& row = (*A_)[inverseOrdering[iter.index()]];
995  for(auto col = row.begin(), cend = row.end(); col != cend; ++col)
996  {
997  iter.insert(ordering_[col.index()]);
998  }
999  }
1000  // Copy values
1001  for(auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
1002  {
1003  auto& newRow = newA[ordering_[iter.index()]];
1004  for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
1005  {
1006  newRow[ordering_[col.index()]] = *col;
1007  }
1008  }
1009  }
1010 
1011  switch ( milu_ )
1012  {
1013  case MILU_VARIANT::MILU_1:
1014  detail::milu0_decomposition ( *ILU);
1015  break;
1016  case MILU_VARIANT::MILU_2:
1017  detail::milu0_decomposition ( *ILU, detail::IdentityFunctor(),
1018  detail::SignFunctor() );
1019  break;
1020  case MILU_VARIANT::MILU_3:
1021  detail::milu0_decomposition ( *ILU, detail::AbsFunctor(),
1022  detail::SignFunctor() );
1023  break;
1024  case MILU_VARIANT::MILU_4:
1025  detail::milu0_decomposition ( *ILU, detail::IdentityFunctor(),
1026  detail::IsPositiveFunctor() );
1027  break;
1028  default:
1029  if (interiorSize_ == A_->N())
1030 #if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
1031  bilu0_decomposition( *ILU );
1032 #else
1033  Dune::ILU::blockILU0Decomposition( *ILU );
1034 #endif
1035  else
1036  detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
1037  break;
1038  }
1039  }
1040  else {
1041  // create ILU-n decomposition
1042  ILU.reset( new Matrix( A_->N(), A_->M(), Matrix::row_wise) );
1043  std::unique_ptr<detail::Reorderer> reorderer, inverseReorderer;
1044  if ( ordering_.empty() )
1045  {
1046  reorderer.reset(new detail::NoReorderer());
1047  inverseReorderer.reset(new detail::NoReorderer());
1048  }
1049  else
1050  {
1051  reorderer.reset(new detail::RealReorderer(ordering_));
1052  inverseReorderer.reset(new detail::RealReorderer(inverseOrdering));
1053  }
1054 
1055  milun_decomposition( *A_, iluIteration_, milu_, *ILU, *reorderer, *inverseReorderer );
1056  }
1057  }
1058  catch (const Dune::MatrixBlockError& error)
1059  {
1060  message = error.what();
1061  std::cerr<<"Exception occurred on process " << rank << " during " <<
1062  "setup of ILU0 preconditioner with message: " <<
1063  message<<std::endl;
1064  ilu_setup_successful = 0;
1065  }
1066 
1067  // Check whether there was a problem on some process
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 )
1071  {
1072  throw Dune::MatrixBlockError();
1073  }
1074 
1075  // store ILU in simple CRS format
1076  detail::convertToCRS( *ILU, lower_, upper_, inv_ );
1077  }
1078 
1079 protected:
1081  Range& reorderD(const Range& d)
1082  {
1083  if ( ordering_.empty())
1084  {
1085  // As d is non-const in the apply method of the
1086  // solver casting away constness in this particular
1087  // setting is not undefined. It is ugly though but due
1088  // to the preconditioner interface of dune-istl.
1089  return const_cast<Range&>(d);
1090  }
1091  else
1092  {
1093  reorderedD_.resize(d.size());
1094  std::size_t i = 0;
1095  for(auto index: ordering_)
1096  {
1097  reorderedD_[index]=d[i++];
1098  }
1099  return reorderedD_;
1100  }
1101  }
1102 
1104  Domain& reorderV(Domain& v)
1105  {
1106  if ( ordering_.empty())
1107  {
1108  return v;
1109  }
1110  else
1111  {
1112  reorderedV_.resize(v.size());
1113  std::size_t i = 0;
1114  for(auto index: ordering_)
1115  {
1116  reorderedV_[index]=v[i++];
1117  }
1118  return reorderedV_;
1119  }
1120  }
1121 
1122  void reorderBack(const Range& reorderedV, Range& v)
1123  {
1124  if ( !ordering_.empty() )
1125  {
1126  std::size_t i = 0;
1127  for(auto index: ordering_)
1128  {
1129  v[i++] = reorderedV[index];
1130  }
1131  }
1132  }
1133 protected:
1136  CRS upper_;
1137  std::vector< block_type > inv_;
1139  std::vector< std::size_t > ordering_;
1143  Domain reorderedV_;
1144 
1145  const ParallelInfo* comm_;
1148  const bool relaxation_;
1149  size_type interiorSize_;
1150  const Matrix* A_;
1151  int iluIteration_;
1152  MILU_VARIANT milu_;
1153  bool redBlack_;
1154  bool reorderSphere_;
1155 };
1156 
1157 } // end namespace Opm
1158 #endif
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