23#include <dune/common/version.hh>
25#include <dune/istl/ilu.hh>
26#include <dune/istl/owneroverlapcopy.hh>
28#include <opm/common/ErrorMacros.hpp>
29#include <opm/common/TimingMacros.hpp>
44 using rowiterator =
typename M::RowIterator;
45 using coliterator =
typename M::ColIterator;
46 using block =
typename M::block_type;
49 for (rowiterator i = A.begin(); i.index() < interiorSize; ++i)
52 coliterator endij=(*i).end();
56 for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
59 coliterator jj = A[ij.index()].find(ij.index());
62 (*ij).rightmultiply(*jj);
65 coliterator endjk=A[ij.index()].end();
66 coliterator jk=jj; ++jk;
67 coliterator ik=ij; ++ik;
68 while (ik!=endij && jk!=endjk)
69 if (ik.index()==jk.index())
78 if (ik.index()<jk.index())
86 if (ij.index()!=i.index())
87 DUNE_THROW(Dune::ISTLError,
"diagonal entry missing");
91 catch (Dune::FMatrixError & e) {
92 DUNE_THROW(Dune::ISTLError,
"ILU failed to invert matrix block");
98template<
class M,
class CRS,
class InvVector>
99void convertToCRS(
const M& A, CRS& lower, CRS& upper, InvVector& inv)
109 using size_type =
typename M :: size_type;
114 lower.resize( A.N() );
115 upper.resize( A.N() );
119 size_type numLower = 0;
120 size_type numUpper = 0;
121 const auto endi = A.end();
122 for (
auto i = A.begin(); i != endi; ++i) {
123 const size_type iIndex = i.index();
124 size_type numLowerRow = 0;
125 for (
auto j = (*i).begin(); j.index() < iIndex; ++j) {
128 numLower += numLowerRow;
129 numUpper += (*i).size() - numLowerRow - 1;
131 assert(numLower + numUpper + A.N() == A.nonzeroes());
133 lower.reserveAdditional( numLower );
137 size_type colcount = 0;
138 lower.rows_[ 0 ] = colcount;
139 for (
auto i=A.begin(); i!=endi; ++i, ++row)
141 const size_type iIndex = i.index();
144 for (
auto j=(*i).begin(); j.index() < iIndex; ++j )
146 lower.push_back( (*j), j.index() );
149 lower.rows_[ iIndex+1 ] = colcount;
152 assert(colcount == numLower);
154 const auto rendi = A.beforeBegin();
157 upper.rows_[ 0 ] = colcount ;
159 upper.reserveAdditional( numUpper );
163 for (
auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
165 const size_type iIndex = i.index();
169 for (
auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
171 const size_type jIndex = j.index();
172 if( j.index() == iIndex )
177 else if ( j.index() >= i.index() )
179 upper.push_back( (*j), jIndex );
183 upper.rows_[ row+1 ] = colcount;
185 assert(colcount == numUpper);
189size_t set_interiorSize( [[maybe_unused]]
size_t N,
size_t interiorSize, [[maybe_unused]]
const PI& comm)
196size_t set_interiorSize(
size_t N,
size_t interiorSize,
const Dune::OwnerOverlapCopyCommunication<int,int>& comm)
200 auto indexSet = comm.indexSet();
203 for (
auto idx = indexSet.begin(); idx!=indexSet.end(); ++idx)
205 if (idx->local().attribute()==1)
207 auto loc = idx->local().local();
220template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
221Dune::SolverCategory::Category
224 return std::is_same_v<ParallelInfoT, Dune::Amg::SequentialInformation> ?
225 Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
228template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
237 comm_(nullptr), w_(w),
238 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
239 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
240 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
248template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
251 const ParallelInfo& comm,
const int n,
const field_type w,
258 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
259 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
260 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
268template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
276template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
286 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
287 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
288 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
296template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
299 const ParallelInfo& comm,
307 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
308 interiorSize_(interiorSize),
309 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
310 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
317template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
319apply (Domain& v,
const Range& d)
321 OPM_TIMEBLOCK(apply);
322 Range& md = reorderD(d);
323 Domain& mv = reorderV(v);
326 using dblock =
typename Range ::block_type;
327 using vblock =
typename Domain::block_type;
331 size_type upperLoopStart = iEnd - interiorSize_;
333 if (iEnd != upper_.rows())
335 OPM_THROW(std::logic_error,
"ILU: number of lower and upper rows must be the same");
339 for (
size_type i = 0; i < lowerLoopEnd; ++i)
341 dblock rhs( md[ i ] );
342 const size_type rowI = lower_.rows_[ i ];
343 const size_type rowINext = lower_.rows_[ i+1 ];
345 for (
size_type col = rowI; col < rowINext; ++col)
347 lower_.values_[ col ].mmv( mv[ lower_.cols_[ col ] ], rhs );
353 for (
size_type i = upperLoopStart; i < iEnd; ++i)
355 vblock& vBlock = mv[ lastRow - i ];
356 vblock rhs ( vBlock );
357 const size_type rowI = upper_.rows_[ i ];
358 const size_type rowINext = upper_.rows_[ i+1 ];
360 for (
size_type col = rowI; col < rowINext; ++col)
362 upper_.values_[ col ].mmv( mv[ upper_.cols_[ col ] ], rhs );
366 inv_[ i ].mv( rhs, vBlock);
369 copyOwnerToAll( mv );
377template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
383 comm_->copyOwnerToAll(v, v);
387template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
391 OPM_TIMEBLOCK(update);
395 if (comm_ && comm_->communicator().size() <= 0)
399 OPM_THROW(std::logic_error,
"Expected a matrix with zero rows for an invalid communicator.");
408 int ilu_setup_successful = 1;
410 const int rank = comm_ ? comm_->communicator().rank() : 0;
414 using Graph = Dune::Amg::MatrixGraph<const Matrix>;
417 const auto& colors = std::get<0>(colorsTuple);
418 const auto& verticesPerColor = std::get<2>(colorsTuple);
419 auto noColors = std::get<1>(colorsTuple);
420 if ( reorderSphere_ )
432 std::vector<std::size_t> inverseOrdering(ordering_.size());
433 std::size_t index = 0;
434 for (
const auto newIndex : ordering_)
436 inverseOrdering[newIndex] = index++;
441 OPM_TIMEBLOCK(iluDecomposition);
442 if (iluIteration_ == 0) {
449 if (ordering_.empty())
452 OPM_TIMEBLOCK(iluDecompositionMakeMatrix);
456 for (std::size_t row = 0; row < A_->N(); ++row) {
457 const auto& Arow = (*A_)[row];
458 auto& ILUrow = (*ILU_)[row];
459 auto Ait = Arow.begin();
460 auto Iit = ILUrow.begin();
461 for (; Ait != Arow.end(); ++Ait, ++Iit) {
467 ILU_ = std::make_unique<Matrix>(*A_);
472 ILU_ = std::make_unique<Matrix>(A_->N(), A_->M(),
473 A_->nonzeroes(), Matrix::row_wise);
476 for (
auto iter = newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
478 const auto& row = (*A_)[inverseOrdering[iter.index()]];
479 for (
auto col = row.begin(), cend = row.end(); col != cend; ++col)
481 iter.insert(ordering_[col.index()]);
485 for (
auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
487 auto& newRow = newA[ordering_[iter.index()]];
488 for (
auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
490 newRow[ordering_[col.index()]] = *col;
502 detail::signFunctor<typename Matrix::field_type> );
506 detail::signFunctor<typename Matrix::field_type> );
510 detail::isPositiveFunctor<typename Matrix::field_type> );
513 if (interiorSize_ == A_->N())
514#if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
515 bilu0_decomposition( *ILU_ );
517 Dune::ILU::blockILU0Decomposition( *ILU_ );
526 ILU_ = std::make_unique<Matrix>(A_->N(), A_->M(), Matrix::row_wise);
527 std::unique_ptr<detail::Reorderer> reorderer, inverseReorderer;
528 if (ordering_.empty())
542 catch (
const Dune::MatrixBlockError& error)
544 message = error.what();
545 std::cerr <<
"Exception occurred on process " << rank <<
" during " <<
546 "setup of ILU0 preconditioner with message: "
547 << message<<std::endl;
548 ilu_setup_successful = 0;
552 const bool parallel_failure = comm_ && comm_->communicator().min(ilu_setup_successful) == 0;
553 const bool local_failure = ilu_setup_successful == 0;
554 if (local_failure || parallel_failure)
556 throw Dune::MatrixBlockError();
563template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
567 if (ordering_.empty())
573 return const_cast<Range&
>(d);
577 reorderedD_.resize(d.size());
579 for (
const auto index : ordering_)
581 reorderedD_[index] = d[i++];
587template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
591 if (ordering_.empty())
597 reorderedV_.resize(v.size());
599 for (
const auto index : ordering_)
601 reorderedV_[index] = v[i++];
607template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
611 if (!ordering_.empty())
614 for (
const auto index : ordering_)
616 v[i++] = reorderedV[index];
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:131
ParallelOverlappingILU0(const Matrix &A, const int n, const field_type w, MILU_VARIANT milu, bool redblack=false, bool reorder_sphere=true)
Constructor.
Definition: ParallelOverlappingILU0_impl.hpp:230
Domain & reorderV(Domain &v)
Reorder V if needed and return a reference to it.
Definition: ParallelOverlappingILU0_impl.hpp:589
size_type interiorSize_
Definition: ParallelOverlappingILU0.hpp:350
void reorderBack(const Range &reorderedV, Range &v)
Definition: ParallelOverlappingILU0_impl.hpp:609
Range & reorderD(const Range &d)
Reorder D if needed and return a reference to it.
Definition: ParallelOverlappingILU0_impl.hpp:565
void copyOwnerToAll(V &v) const
Definition: ParallelOverlappingILU0_impl.hpp:380
Dune::SolverCategory::Category category() const override
Definition: ParallelOverlappingILU0_impl.hpp:222
void update() override
Definition: ParallelOverlappingILU0_impl.hpp:389
void apply(Domain &v, const Range &d) override
Apply the preconditoner.
Definition: ParallelOverlappingILU0_impl.hpp:319
typename matrix_type::size_type size_type
Definition: ParallelOverlappingILU0.hpp:145
typename Domain::field_type field_type
The field type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:142
void milu0_decomposition(M &A, FieldFunct< M > absFunctor=signFunctor< typename M::field_type >, FieldFunct< M > signFunctor=oneFunctor< typename M::field_type >, std::vector< typename M::block_type > *diagonal=nullptr)
void convertToCRS(const M &A, CRS &lower, CRS &upper, InvVector &inv)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ParallelOverlappingILU0_impl.hpp:99
void milun_decomposition(const M &A, int n, MILU_VARIANT milu, M &ILU, Reorderer &ordering, Reorderer &inverseOrdering)
size_t set_interiorSize(size_t N, size_t interiorSize, const PI &comm)
Definition: ParallelOverlappingILU0_impl.hpp:189
void ghost_last_bilu0_decomposition(M &A, std::size_t interiorSize)
Compute Blocked ILU0 decomposition, when we know junk ghost rows are located at the end of A.
Definition: ParallelOverlappingILU0_impl.hpp:41
Definition: blackoilboundaryratevector.hh:37
MILU_VARIANT
Definition: MILU.hpp:34
@ MILU_1
sum(dropped entries)
@ MILU_2
sum(dropped entries)
@ MILU_3
sum(|dropped entries|)
@ MILU_4
sum(dropped entries)
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:169
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:189
std::tuple< std::vector< int >, int, std::vector< std::size_t > > colorVerticesWelshPowell(const Graph &graph)
Color the vertices of graph.
Definition: GraphColoring.hpp:113