My Project
ParallelOverlappingILU0_impl.hpp
1/*
2 Copyright 2015, 2022 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
21#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
22
23#include <dune/istl/ilu.hh>
24#include <dune/istl/owneroverlapcopy.hh>
25
26#include <opm/common/ErrorMacros.hpp>
27
28#include <opm/simulators/linalg/GraphColoring.hpp>
29#include <opm/simulators/linalg/matrixblock.hh>
30
31namespace Opm
32{
33namespace detail
34{
35
37template<class M>
38void ghost_last_bilu0_decomposition (M& A, size_t interiorSize)
39{
40 // iterator types
41 using rowiterator = typename M::RowIterator;
42 using coliterator = typename M::ColIterator;
43 using block = typename M::block_type;
44
45 // implement left looking variant with stored inverse
46 for (rowiterator i = A.begin(); i.index() < interiorSize; ++i)
47 {
48 // coliterator is diagonal after the following loop
49 coliterator endij=(*i).end(); // end of row i
50 coliterator ij;
51
52 // eliminate entries left of diagonal; store L factor
53 for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
54 {
55 // find A_jj which eliminates A_ij
56 coliterator jj = A[ij.index()].find(ij.index());
57
58 // compute L_ij = A_jj^-1 * A_ij
59 (*ij).rightmultiply(*jj);
60
61 // modify row
62 coliterator endjk=A[ij.index()].end(); // end of row j
63 coliterator jk=jj; ++jk;
64 coliterator ik=ij; ++ik;
65 while (ik!=endij && jk!=endjk)
66 if (ik.index()==jk.index())
67 {
68 block B(*jk);
69 B.leftmultiply(*ij);
70 *ik -= B;
71 ++ik; ++jk;
72 }
73 else
74 {
75 if (ik.index()<jk.index())
76 ++ik;
77 else
78 ++jk;
79 }
80 }
81
82 // invert pivot and store it in A
83 if (ij.index()!=i.index())
84 DUNE_THROW(Dune::ISTLError,"diagonal entry missing");
85 try {
86 (*ij).invert(); // compute inverse of diagonal block
87 }
88 catch (Dune::FMatrixError & e) {
89 DUNE_THROW(Dune::ISTLError,"ILU failed to invert matrix block");
90 }
91 }
92}
93
95template<class M, class CRS, class InvVector>
96void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv)
97{
98 // No need to do anything for 0 rows. Return to prevent indexing a
99 // a zero sized array.
100 if ( A.N() == 0 )
101 {
102 return;
103 }
104
105 using size_type = typename M :: size_type;
106
107 lower.clear();
108 upper.clear();
109 inv.clear();
110 lower.resize( A.N() );
111 upper.resize( A.N() );
112 inv.resize( A.N() );
113
114 // Count the lower and upper matrix entries.
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) {
122 ++numLowerRow;
123 }
124 numLower += numLowerRow;
125 numUpper += (*i).size() - numLowerRow - 1;
126 }
127 assert(numLower + numUpper + A.N() == A.nonzeroes());
128
129 lower.reserveAdditional( numLower );
130
131 // implement left looking variant with stored inverse
132 size_type row = 0;
133 size_type colcount = 0;
134 lower.rows_[ 0 ] = colcount;
135 for (auto i=A.begin(); i!=endi; ++i, ++row)
136 {
137 const size_type iIndex = i.index();
138
139 // eliminate entries left of diagonal; store L factor
140 for (auto j=(*i).begin(); j.index() < iIndex; ++j )
141 {
142 lower.push_back( (*j), j.index() );
143 ++colcount;
144 }
145 lower.rows_[ iIndex+1 ] = colcount;
146 }
147
148 assert(colcount == numLower);
149
150 const auto rendi = A.beforeBegin();
151 row = 0;
152 colcount = 0;
153 upper.rows_[ 0 ] = colcount ;
154
155 upper.reserveAdditional( numUpper );
156
157 // NOTE: upper and inv store entries in reverse order, reverse here
158 // relative to ILU
159 for (auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
160 {
161 const size_type iIndex = i.index();
162
163 // store in reverse row order
164 // eliminate entries left of diagonal; store L factor
165 for (auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
166 {
167 const size_type jIndex = j.index();
168 if( j.index() == iIndex )
169 {
170 inv[ row ] = (*j);
171 break;
172 }
173 else if ( j.index() >= i.index() )
174 {
175 upper.push_back( (*j), jIndex );
176 ++colcount ;
177 }
178 }
179 upper.rows_[ row+1 ] = colcount;
180 }
181 assert(colcount == numUpper);
182}
183
184} // end namespace detail
185
186
187template<class Matrix, class Domain, class Range, class ParallelInfoT>
188Dune::SolverCategory::Category
189ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfoT>::category() const
190{
191 return std::is_same_v<ParallelInfoT, Dune::Amg::SequentialInformation> ?
192 Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
193}
194
195template<class Matrix, class Domain, class Range, class ParallelInfoT>
197ParallelOverlappingILU0(const Matrix& A,
198 const int n, const field_type w,
199 MILU_VARIANT milu, bool redblack,
200 bool reorder_sphere)
201 : lower_(),
202 upper_(),
203 inv_(),
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)
208{
209 interiorSize_ = A.N();
210 // BlockMatrix is a Subclass of FieldMatrix that just adds
211 // methods. Therefore this cast should be safe.
212 update();
213}
214
215template<class Matrix, class Domain, class Range, class ParallelInfoT>
217ParallelOverlappingILU0(const Matrix& A,
218 const ParallelInfo& comm, const int n, const field_type w,
219 MILU_VARIANT milu, bool redblack,
220 bool reorder_sphere)
221 : lower_(),
222 upper_(),
223 inv_(),
224 comm_(&comm), w_(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)
228{
229 interiorSize_ = A.N();
230 // BlockMatrix is a Subclass of FieldMatrix that just adds
231 // methods. Therefore this cast should be safe.
232 update();
233}
234
235template<class Matrix, class Domain, class Range, class ParallelInfoT>
237ParallelOverlappingILU0(const Matrix& A,
238 const field_type w, MILU_VARIANT milu, bool redblack,
239 bool reorder_sphere)
240 : ParallelOverlappingILU0( A, 0, w, milu, redblack, reorder_sphere )
241{}
242
243template<class Matrix, class Domain, class Range, class ParallelInfoT>
245ParallelOverlappingILU0(const Matrix& A,
246 const ParallelInfo& comm, const field_type w,
247 MILU_VARIANT milu, bool redblack,
248 bool reorder_sphere)
249 : lower_(),
250 upper_(),
251 inv_(),
252 comm_(&comm), w_(w),
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)
256{
257 interiorSize_ = A.N();
258 // BlockMatrix is a Subclass of FieldMatrix that just adds
259 // methods. Therefore this cast should be safe.
260 update();
261}
262
263template<class Matrix, class Domain, class Range, class ParallelInfoT>
265ParallelOverlappingILU0(const Matrix& A,
266 const ParallelInfo& comm,
267 const field_type w, MILU_VARIANT milu,
268 size_type interiorSize, bool redblack,
269 bool reorder_sphere)
270 : lower_(),
271 upper_(),
272 inv_(),
273 comm_(&comm), w_(w),
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)
278{
279 // BlockMatrix is a Subclass of FieldMatrix that just adds
280 // methods. Therefore this cast should be safe.
281 update( );
282}
283
284template<class Matrix, class Domain, class Range, class ParallelInfoT>
286apply (Domain& v, const Range& d)
287{
288 Range& md = reorderD(d);
289 Domain& mv = reorderV(v);
290
291 // iterator types
292 using dblock = typename Range ::block_type;
293 using vblock = typename Domain::block_type;
294
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())
300 {
301 OPM_THROW(std::logic_error,"ILU: number of lower and upper rows must be the same");
302 }
303
304 // lower triangular solve
305 for (size_type i = 0; i < lowerLoopEnd; ++i)
306 {
307 dblock rhs( md[ i ] );
308 const size_type rowI = lower_.rows_[ i ];
309 const size_type rowINext = lower_.rows_[ i+1 ];
310
311 for (size_type col = rowI; col < rowINext; ++col)
312 {
313 lower_.values_[ col ].mmv( mv[ lower_.cols_[ col ] ], rhs );
314 }
315
316 mv[ i ] = rhs; // Lii = I
317 }
318
319 for (size_type i = upperLoopStart; i < iEnd; ++i)
320 {
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 ];
325
326 for (size_type col = rowI; col < rowINext; ++col)
327 {
328 upper_.values_[ col ].mmv( mv[ upper_.cols_[ col ] ], rhs );
329 }
330
331 // apply inverse and store result
332 inv_[ i ].mv( rhs, vBlock);
333 }
334
335 copyOwnerToAll( mv );
336
337 if( relaxation_ ) {
338 mv *= w_;
339 }
340 reorderBack(mv, v);
341}
342
343template<class Matrix, class Domain, class Range, class ParallelInfoT>
344template<class V>
346copyOwnerToAll(V& v) const
347{
348 if( comm_ ) {
349 comm_->copyOwnerToAll(v, v);
350 }
351}
352
353template<class Matrix, class Domain, class Range, class ParallelInfoT>
354void ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfoT>::
355update()
356{
357 // (For older DUNE versions the communicator might be
358 // invalid if redistribution in AMG happened on the coarset level.
359 // Therefore we check for nonzero size
360 if (comm_ && comm_->communicator().size() <= 0)
361 {
362 if (A_->N() > 0)
363 {
364 OPM_THROW(std::logic_error, "Expected a matrix with zero rows for an invalid communicator.");
365 }
366 else
367 {
368 // simply set the communicator to null
369 comm_ = nullptr;
370 }
371 }
372
373 int ilu_setup_successful = 1;
374 std::string message;
375 const int rank = comm_ ? comm_->communicator().rank() : 0;
376
377 std::unique_ptr<Matrix> ILU;
378
379 if (redBlack_)
380 {
381 using Graph = Dune::Amg::MatrixGraph<const Matrix>;
382 Graph graph(*A_);
383 auto colorsTuple = colorVerticesWelshPowell(graph);
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_ )
388 {
389 ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor,
390 graph, 0);
391 }
392 else
393 {
394 ordering_ = reorderVerticesPreserving(colors, noColors, verticesPerColor,
395 graph);
396 }
397 }
398
399 std::vector<std::size_t> inverseOrdering(ordering_.size());
400 std::size_t index = 0;
401 for (const auto newIndex : ordering_)
402 {
403 inverseOrdering[newIndex] = index++;
404 }
405
406 try
407 {
408 if (iluIteration_ == 0) {
409 // create ILU-0 decomposition
410 if (ordering_.empty())
411 {
412 ILU = std::make_unique<Matrix>(*A_);
413 }
414 else
415 {
416 ILU = std::make_unique<Matrix>(A_->N(), A_->M(),
417 A_->nonzeroes(), Matrix::row_wise);
418 auto& newA = *ILU;
419 // Create sparsity pattern
420 for (auto iter = newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
421 {
422 const auto& row = (*A_)[inverseOrdering[iter.index()]];
423 for (auto col = row.begin(), cend = row.end(); col != cend; ++col)
424 {
425 iter.insert(ordering_[col.index()]);
426 }
427 }
428 // Copy values
429 for (auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
430 {
431 auto& newRow = newA[ordering_[iter.index()]];
432 for (auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
433 {
434 newRow[ordering_[col.index()]] = *col;
435 }
436 }
437 }
438
439 switch (milu_)
440 {
442 detail::milu0_decomposition ( *ILU);
443 break;
445 detail::milu0_decomposition ( *ILU, detail::identityFunctor<typename Matrix::field_type>,
446 detail::signFunctor<typename Matrix::field_type> );
447 break;
449 detail::milu0_decomposition ( *ILU, detail::absFunctor<typename Matrix::field_type>,
450 detail::signFunctor<typename Matrix::field_type> );
451 break;
453 detail::milu0_decomposition ( *ILU, detail::identityFunctor<typename Matrix::field_type>,
454 detail::isPositiveFunctor<typename Matrix::field_type> );
455 break;
456 default:
457 if (interiorSize_ == A_->N())
458#if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
459 bilu0_decomposition( *ILU );
460#else
461 Dune::ILU::blockILU0Decomposition( *ILU );
462#endif
463 else
464 detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
465 break;
466 }
467 }
468 else {
469 // create ILU-n decomposition
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())
473 {
474 reorderer.reset(new detail::NoReorderer());
475 inverseReorderer.reset(new detail::NoReorderer());
476 }
477 else
478 {
479 reorderer.reset(new detail::RealReorderer(ordering_));
480 inverseReorderer.reset(new detail::RealReorderer(inverseOrdering));
481 }
482
483 milun_decomposition( *A_, iluIteration_, milu_, *ILU, *reorderer, *inverseReorderer );
484 }
485 }
486 catch (const Dune::MatrixBlockError& error)
487 {
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;
493 }
494
495 // Check whether there was a problem on some process
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)
499 {
500 throw Dune::MatrixBlockError();
501 }
502
503 // store ILU in simple CRS format
504 detail::convertToCRS(*ILU, lower_, upper_, inv_);
505}
506
507template<class Matrix, class Domain, class Range, class ParallelInfoT>
509reorderD(const Range& d)
510{
511 if (ordering_.empty())
512 {
513 // As d is non-const in the apply method of the
514 // solver casting away constness in this particular
515 // setting is not undefined. It is ugly though but due
516 // to the preconditioner interface of dune-istl.
517 return const_cast<Range&>(d);
518 }
519 else
520 {
521 reorderedD_.resize(d.size());
522 std::size_t i = 0;
523 for (const auto index : ordering_)
524 {
525 reorderedD_[index] = d[i++];
526 }
527 return reorderedD_;
528 }
529}
530
531template<class Matrix, class Domain, class Range, class ParallelInfoT>
533reorderV(Domain& v)
534{
535 if (ordering_.empty())
536 {
537 return v;
538 }
539 else
540 {
541 reorderedV_.resize(v.size());
542 std::size_t i = 0;
543 for (const auto index : ordering_)
544 {
545 reorderedV_[index] = v[i++];
546 }
547 return reorderedV_;
548 }
549}
550
551template<class Matrix, class Domain, class Range, class ParallelInfoT>
553reorderBack(const Range& reorderedV, Range& v)
554{
555 if (!ordering_.empty())
556 {
557 std::size_t i = 0;
558 for (const auto index : ordering_)
559 {
560 v[i++] = reorderedV[index];
561 }
562 }
563}
564
565} // end namespace Opm
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