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>
32#include <opm/simulators/linalg/matrixblock.hh>
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);
191template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
192Dune::SolverCategory::Category
195 return std::is_same_v<ParallelInfoT, Dune::Amg::SequentialInformation> ?
196 Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
199template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
208 comm_(nullptr), w_(w),
209 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
210 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
211 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
219template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
222 const ParallelInfo& comm,
const int n,
const field_type w,
229 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
230 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
231 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
239template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
247template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
257 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
258 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
259 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
267template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
270 const ParallelInfo& comm,
278 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
279 interiorSize_(interiorSize),
280 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
281 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
288template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
290apply (Domain& v,
const Range& d)
292 OPM_TIMEBLOCK(apply);
293 Range& md = reorderD(d);
294 Domain& mv = reorderV(v);
297 using dblock =
typename Range ::block_type;
298 using vblock =
typename Domain::block_type;
302 size_type upperLoopStart = iEnd - interiorSize_;
304 if (iEnd != upper_.rows())
306 OPM_THROW(std::logic_error,
"ILU: number of lower and upper rows must be the same");
310 for (
size_type i = 0; i < lowerLoopEnd; ++i)
312 dblock rhs( md[ i ] );
313 const size_type rowI = lower_.rows_[ i ];
314 const size_type rowINext = lower_.rows_[ i+1 ];
316 for (
size_type col = rowI; col < rowINext; ++col)
318 lower_.values_[ col ].mmv( mv[ lower_.cols_[ col ] ], rhs );
324 for (
size_type i = upperLoopStart; i < iEnd; ++i)
326 vblock& vBlock = mv[ lastRow - i ];
327 vblock rhs ( vBlock );
328 const size_type rowI = upper_.rows_[ i ];
329 const size_type rowINext = upper_.rows_[ i+1 ];
331 for (
size_type col = rowI; col < rowINext; ++col)
333 upper_.values_[ col ].mmv( mv[ upper_.cols_[ col ] ], rhs );
337 inv_[ i ].mv( rhs, vBlock);
340 copyOwnerToAll( mv );
348template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
354 comm_->copyOwnerToAll(v, v);
358template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
362 OPM_TIMEBLOCK(update);
366 if (comm_ && comm_->communicator().size() <= 0)
370 OPM_THROW(std::logic_error,
"Expected a matrix with zero rows for an invalid communicator.");
379 int ilu_setup_successful = 1;
381 const int rank = comm_ ? comm_->communicator().rank() : 0;
385 using Graph = Dune::Amg::MatrixGraph<const Matrix>;
388 const auto& colors = std::get<0>(colorsTuple);
389 const auto& verticesPerColor = std::get<2>(colorsTuple);
390 auto noColors = std::get<1>(colorsTuple);
391 if ( reorderSphere_ )
403 std::vector<std::size_t> inverseOrdering(ordering_.size());
404 std::size_t index = 0;
405 for (
const auto newIndex : ordering_)
407 inverseOrdering[newIndex] = index++;
412 OPM_TIMEBLOCK(iluDecomposition);
413 if (iluIteration_ == 0) {
415 if (ordering_.empty())
418 OPM_TIMEBLOCK(iluDecompositionMakeMatrix);
422 for (std::size_t row = 0; row < A_->N(); ++row) {
423 const auto& Arow = (*A_)[row];
424 auto& ILUrow = (*ILU_)[row];
425 auto Ait = Arow.begin();
426 auto Iit = ILUrow.begin();
427 for (; Ait != Arow.end(); ++Ait, ++Iit) {
433 ILU_ = std::make_unique<Matrix>(*A_);
438 ILU_ = std::make_unique<Matrix>(A_->N(), A_->M(),
439 A_->nonzeroes(), Matrix::row_wise);
442 for (
auto iter = newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
444 const auto& row = (*A_)[inverseOrdering[iter.index()]];
445 for (
auto col = row.begin(), cend = row.end(); col != cend; ++col)
447 iter.insert(ordering_[col.index()]);
451 for (
auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
453 auto& newRow = newA[ordering_[iter.index()]];
454 for (
auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
456 newRow[ordering_[col.index()]] = *col;
468 detail::signFunctor<typename Matrix::field_type> );
472 detail::signFunctor<typename Matrix::field_type> );
476 detail::isPositiveFunctor<typename Matrix::field_type> );
479 if (interiorSize_ == A_->N())
480#if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
481 bilu0_decomposition( *ILU_ );
483 Dune::ILU::blockILU0Decomposition( *ILU_ );
492 ILU_ = std::make_unique<Matrix>(A_->N(), A_->M(), Matrix::row_wise);
493 std::unique_ptr<detail::Reorderer> reorderer, inverseReorderer;
494 if (ordering_.empty())
508 catch (
const Dune::MatrixBlockError& error)
510 message = error.what();
511 std::cerr <<
"Exception occurred on process " << rank <<
" during " <<
512 "setup of ILU0 preconditioner with message: "
513 << message<<std::endl;
514 ilu_setup_successful = 0;
518 const bool parallel_failure = comm_ && comm_->communicator().min(ilu_setup_successful) == 0;
519 const bool local_failure = ilu_setup_successful == 0;
520 if (local_failure || parallel_failure)
522 throw Dune::MatrixBlockError();
529template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
533 if (ordering_.empty())
539 return const_cast<Range&
>(d);
543 reorderedD_.resize(d.size());
545 for (
const auto index : ordering_)
547 reorderedD_[index] = d[i++];
553template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
557 if (ordering_.empty())
563 reorderedV_.resize(v.size());
565 for (
const auto index : ordering_)
567 reorderedV_[index] = v[i++];
573template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
577 if (!ordering_.empty())
580 for (
const auto index : ordering_)
582 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:201
Domain & reorderV(Domain &v)
Reorder V if needed and return a reference to it.
Definition: ParallelOverlappingILU0_impl.hpp:555
size_type interiorSize_
Definition: ParallelOverlappingILU0.hpp:346
void reorderBack(const Range &reorderedV, Range &v)
Definition: ParallelOverlappingILU0_impl.hpp:575
Range & reorderD(const Range &d)
Reorder D if needed and return a reference to it.
Definition: ParallelOverlappingILU0_impl.hpp:531
void copyOwnerToAll(V &v) const
Definition: ParallelOverlappingILU0_impl.hpp:351
Dune::SolverCategory::Category category() const override
Definition: ParallelOverlappingILU0_impl.hpp:193
void update() override
Definition: ParallelOverlappingILU0_impl.hpp:360
void apply(Domain &v, const Range &d) override
Apply the preconditoner.
Definition: ParallelOverlappingILU0_impl.hpp:290
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)
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: BlackoilPhases.hpp:27
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