21 #include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
23 #include <dune/istl/ilu.hh>
24 #include <dune/istl/owneroverlapcopy.hh>
26 #include <opm/common/ErrorMacros.hpp>
28 #include <opm/simulators/linalg/GraphColoring.hpp>
29 #include <opm/simulators/linalg/matrixblock.hh>
38 void ghost_last_bilu0_decomposition (M& A,
size_t interiorSize)
41 using rowiterator =
typename M::RowIterator;
42 using coliterator =
typename M::ColIterator;
43 using block =
typename M::block_type;
46 for (rowiterator i = A.begin(); i.index() < interiorSize; ++i)
49 coliterator endij=(*i).end();
53 for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
56 coliterator jj = A[ij.index()].find(ij.index());
59 (*ij).rightmultiply(*jj);
62 coliterator endjk=A[ij.index()].end();
63 coliterator jk=jj; ++jk;
64 coliterator ik=ij; ++ik;
65 while (ik!=endij && jk!=endjk)
66 if (ik.index()==jk.index())
75 if (ik.index()<jk.index())
83 if (ij.index()!=i.index())
84 DUNE_THROW(Dune::ISTLError,
"diagonal entry missing");
88 catch (Dune::FMatrixError & e) {
89 DUNE_THROW(Dune::ISTLError,
"ILU failed to invert matrix block");
95 template<
class M,
class CRS,
class InvVector>
96 void convertToCRS(
const M& A, CRS& lower, CRS& upper, InvVector& inv)
105 using size_type =
typename M :: size_type;
110 lower.resize( A.N() );
111 upper.resize( A.N() );
115 size_type numLower = 0;
116 size_type numUpper = 0;
117 const auto endi = A.end();
118 for (
auto i = A.begin(); i != endi; ++i) {
119 const size_type iIndex = i.index();
120 size_type numLowerRow = 0;
121 for (
auto j = (*i).begin(); j.index() < iIndex; ++j) {
124 numLower += numLowerRow;
125 numUpper += (*i).size() - numLowerRow - 1;
127 assert(numLower + numUpper + A.N() == A.nonzeroes());
129 lower.reserveAdditional( numLower );
133 size_type colcount = 0;
134 lower.rows_[ 0 ] = colcount;
135 for (
auto i=A.begin(); i!=endi; ++i, ++row)
137 const size_type iIndex = i.index();
140 for (
auto j=(*i).begin(); j.index() < iIndex; ++j )
142 lower.push_back( (*j), j.index() );
145 lower.rows_[ iIndex+1 ] = colcount;
148 assert(colcount == numLower);
150 const auto rendi = A.beforeBegin();
153 upper.rows_[ 0 ] = colcount ;
155 upper.reserveAdditional( numUpper );
159 for (
auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
161 const size_type iIndex = i.index();
165 for (
auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
167 const size_type jIndex = j.index();
168 if( j.index() == iIndex )
173 else if ( j.index() >= i.index() )
175 upper.push_back( (*j), jIndex );
179 upper.rows_[ row+1 ] = colcount;
181 assert(colcount == numUpper);
187 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
188 Dune::SolverCategory::Category
189 ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfoT>::category()
const
191 return std::is_same_v<ParallelInfoT, Dune::Amg::SequentialInformation> ?
192 Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
195 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
204 comm_(nullptr), w_(w),
205 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
206 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
207 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
209 interiorSize_ = A.N();
215 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
218 const ParallelInfo& comm,
const int n,
const field_type w,
225 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
226 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
227 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
229 interiorSize_ = A.N();
235 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
243 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
253 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
254 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
255 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
257 interiorSize_ = A.N();
263 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
266 const ParallelInfo& comm,
268 size_type interiorSize,
bool redblack,
274 relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
275 interiorSize_(interiorSize),
276 A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
277 milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
284 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
286 apply (Domain& v,
const Range& d)
288 Range& md = reorderD(d);
289 Domain& mv = reorderV(v);
292 using dblock =
typename Range ::block_type;
293 using vblock =
typename Domain::block_type;
295 const size_type iEnd = lower_.rows();
296 const size_type lastRow = iEnd - 1;
297 size_type upperLoopStart = iEnd - interiorSize_;
298 size_type lowerLoopEnd = interiorSize_;
299 if (iEnd != upper_.rows())
301 OPM_THROW(std::logic_error,
"ILU: number of lower and upper rows must be the same");
305 for (size_type i = 0; i < lowerLoopEnd; ++i)
307 dblock rhs( md[ i ] );
308 const size_type rowI = lower_.rows_[ i ];
309 const size_type rowINext = lower_.rows_[ i+1 ];
311 for (size_type col = rowI; col < rowINext; ++col)
313 lower_.values_[ col ].mmv( mv[ lower_.cols_[ col ] ], rhs );
319 for (size_type i = upperLoopStart; i < iEnd; ++i)
321 vblock& vBlock = mv[ lastRow - i ];
322 vblock rhs ( vBlock );
323 const size_type rowI = upper_.rows_[ i ];
324 const size_type rowINext = upper_.rows_[ i+1 ];
326 for (size_type col = rowI; col < rowINext; ++col)
328 upper_.values_[ col ].mmv( mv[ upper_.cols_[ col ] ], rhs );
332 inv_[ i ].mv( rhs, vBlock);
335 copyOwnerToAll( mv );
343 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
349 comm_->copyOwnerToAll(v, v);
353 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
354 void ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfoT>::
360 if (comm_ && comm_->communicator().size() <= 0)
364 OPM_THROW(std::logic_error,
"Expected a matrix with zero rows for an invalid communicator.");
373 int ilu_setup_successful = 1;
375 const int rank = comm_ ? comm_->communicator().rank() : 0;
377 std::unique_ptr<Matrix>
ILU;
381 using Graph = Dune::Amg::MatrixGraph<const Matrix>;
384 const auto& colors = std::get<0>(colorsTuple);
385 const auto& verticesPerColor = std::get<2>(colorsTuple);
386 auto noColors = std::get<1>(colorsTuple);
387 if ( reorderSphere_ )
399 std::vector<std::size_t> inverseOrdering(ordering_.size());
400 std::size_t index = 0;
401 for (
const auto newIndex : ordering_)
403 inverseOrdering[newIndex] = index++;
408 if (iluIteration_ == 0) {
410 if (ordering_.empty())
412 ILU = std::make_unique<Matrix>(*A_);
416 ILU = std::make_unique<Matrix>(A_->N(), A_->M(),
417 A_->nonzeroes(), Matrix::row_wise);
420 for (
auto iter = newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
422 const auto& row = (*A_)[inverseOrdering[iter.index()]];
423 for (
auto col = row.begin(), cend = row.end(); col != cend; ++col)
425 iter.insert(ordering_[col.index()]);
429 for (
auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
431 auto& newRow = newA[ordering_[iter.index()]];
432 for (
auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
434 newRow[ordering_[col.index()]] = *col;
442 detail::milu0_decomposition ( *ILU);
445 detail::milu0_decomposition ( *ILU, detail::identityFunctor<typename Matrix::field_type>,
446 detail::signFunctor<typename Matrix::field_type> );
449 detail::milu0_decomposition ( *ILU, detail::absFunctor<typename Matrix::field_type>,
450 detail::signFunctor<typename Matrix::field_type> );
453 detail::milu0_decomposition ( *ILU, detail::identityFunctor<typename Matrix::field_type>,
454 detail::isPositiveFunctor<typename Matrix::field_type> );
457 if (interiorSize_ == A_->N())
458 #if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
459 bilu0_decomposition( *ILU );
461 Dune::ILU::blockILU0Decomposition( *ILU );
464 detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
470 ILU = std::make_unique<Matrix>(A_->N(), A_->M(), Matrix::row_wise);
471 std::unique_ptr<detail::Reorderer> reorderer, inverseReorderer;
472 if (ordering_.empty())
474 reorderer.reset(
new detail::NoReorderer());
475 inverseReorderer.reset(
new detail::NoReorderer());
479 reorderer.reset(
new detail::RealReorderer(ordering_));
480 inverseReorderer.reset(
new detail::RealReorderer(inverseOrdering));
483 milun_decomposition( *A_, iluIteration_, milu_, *ILU, *reorderer, *inverseReorderer );
486 catch (
const Dune::MatrixBlockError& error)
488 message = error.what();
489 std::cerr <<
"Exception occurred on process " << rank <<
" during " <<
490 "setup of ILU0 preconditioner with message: "
491 << message<<std::endl;
492 ilu_setup_successful = 0;
496 const bool parallel_failure = comm_ && comm_->communicator().min(ilu_setup_successful) == 0;
497 const bool local_failure = ilu_setup_successful == 0;
498 if (local_failure || parallel_failure)
500 throw Dune::MatrixBlockError();
504 detail::convertToCRS(*ILU, lower_, upper_, inv_);
507 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
511 if (ordering_.empty())
517 return const_cast<Range&
>(d);
521 reorderedD_.resize(d.size());
523 for (
const auto index : ordering_)
525 reorderedD_[index] = d[i++];
531 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
535 if (ordering_.empty())
541 reorderedV_.resize(v.size());
543 for (
const auto index : ordering_)
545 reorderedV_[index] = v[i++];
551 template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
555 if (!ordering_.empty())
558 for (
const auto index : ordering_)
560 v[i++] = reorderedV[index];
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:145
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:197
Domain & reorderV(Domain &v)
Reorder V if needed and return a reference to it.
Definition: ParallelOverlappingILU0_impl.hpp:533
Range & reorderD(const Range &d)
Reorder D if needed and return a reference to it.
Definition: ParallelOverlappingILU0_impl.hpp:509
void apply(Domain &v, const Range &d) override
Apply the preconditoner.
Definition: ParallelOverlappingILU0_impl.hpp:286
typename Domain::field_type field_type
The field type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:156
This file contains a set of helper functions used by VFPProd / VFPInj.
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)
@ 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