21#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
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>
31#include <opm/simulators/linalg/GraphColoring.hpp>
32#include <opm/simulators/linalg/matrixblock.hh>
41void ghost_last_bilu0_decomposition (M&
A, std::size_t
interiorSize)
46 using block =
typename M::block_type;
56 for (
ij=(*i).begin();
ij.index()<i.index(); ++
ij)
62 (*ij).rightmultiply(*
jj);
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() );
121 const auto endi =
A.end();
122 for (
auto i =
A.begin(); i !=
endi; ++i) {
123 const size_type
iIndex = i.index();
125 for (
auto j = (*i).begin(); j.index() <
iIndex; ++j) {
133 lower.reserveAdditional(
numLower );
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() );
154 const auto rendi =
A.beforeBegin();
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();
177 else if ( j.index() >= i.index() )
179 upper.push_back( (*j),
jIndex );
196size_t set_interiorSize(
size_t N,
size_t interiorSize,
const Dune::OwnerOverlapCopyCommunication<int,int>& comm)
205 if (
idx->local().attribute()==1)
207 auto loc =
idx->local().local();
220template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
221Dune::SolverCategory::Category
222ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfoT>::category()
const
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>
238 relaxation_( std::
abs(
w - 1.0 ) > 1
e-15 ),
242 interiorSize_ =
A.N();
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 ) > 1
e-15 ),
262 interiorSize_ =
A.N();
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 ) > 1
e-15 ),
290 interiorSize_ =
A.N();
296template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
299 const ParallelInfo& comm,
307 relaxation_( std::
abs(
w - 1.0 ) > 1
e-15 ),
317template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
319apply (Domain&
v,
const Range&
d)
322 Range&
md = reorderD(
d);
323 Domain& mv = reorderV(
v);
326 using dblock =
typename Range ::block_type;
327 using vblock =
typename Domain::block_type;
329 const size_type
iEnd = lower_.rows();
333 if (
iEnd != upper_.rows())
335 OPM_THROW(std::logic_error,
"ILU: number of lower and upper rows must be the same");
342 const size_type
rowI = lower_.rows_[ i ];
343 const size_type
rowINext = lower_.rows_[ i+1 ];
347 lower_.values_[
col ].mmv( mv[ lower_.cols_[
col ] ], rhs );
357 const size_type
rowI = upper_.rows_[ i ];
358 const size_type
rowINext = upper_.rows_[ i+1 ];
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>
388void ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfoT>::
395 if (comm_ && comm_->communicator().size() <= 0)
399 OPM_THROW(std::logic_error,
"Expected a matrix with zero rows for an invalid communicator.");
410 const int rank = comm_ ? comm_->communicator().rank() : 0;
414 using Graph = Dune::Amg::MatrixGraph<const Matrix>;
420 if ( reorderSphere_ )
433 std::size_t index = 0;
434 for (
const auto newIndex : ordering_)
442 if (iluIteration_ == 0) {
445 interiorSize_ = detail::set_interiorSize(A_->N(), interiorSize_, *comm_);
449 if (ordering_.empty())
456 for (std::size_t
row = 0;
row < A_->N(); ++
row) {
467 ILU_ = std::make_unique<Matrix>(*A_);
472 ILU_ = std::make_unique<Matrix>(A_->N(), A_->M(),
473 A_->nonzeroes(), Matrix::row_wise);
481 iter.insert(ordering_[
col.index()]);
498 detail::milu0_decomposition ( *ILU_);
501 detail::milu0_decomposition ( *ILU_, detail::identityFunctor<typename Matrix::field_type>,
502 detail::signFunctor<typename Matrix::field_type> );
505 detail::milu0_decomposition ( *ILU_, detail::absFunctor<typename Matrix::field_type>,
506 detail::signFunctor<typename Matrix::field_type> );
509 detail::milu0_decomposition ( *ILU_, detail::identityFunctor<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)
517 Dune::ILU::blockILU0Decomposition( *ILU_ );
520 detail::ghost_last_bilu0_decomposition(*ILU_, interiorSize_);
526 ILU_ = std::make_unique<Matrix>(A_->N(), A_->M(), Matrix::row_wise);
528 if (ordering_.empty())
530 reorderer.reset(
new detail::NoReorderer());
535 reorderer.reset(
new detail::RealReorderer(ordering_));
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;
556 throw Dune::MatrixBlockError();
560 detail::convertToCRS(*ILU_, lower_, upper_, inv_);
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_)
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
Range & reorderD(const Range &d)
Reorder D if needed and return a reference to it.
Definition ParallelOverlappingILU0_impl.hpp:565
void apply(Domain &v, const Range &d) override
Apply the preconditoner.
Definition ParallelOverlappingILU0_impl.hpp:319
typename Domain::field_type field_type
The field type of the preconditioner.
Definition ParallelOverlappingILU0.hpp:142
This file contains a set of helper functions used by VFPProd / VFPInj.
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
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
std::tuple< std::vector< int >, int, std::vector< std::size_t > > colorVerticesWelshPowell(const Graph &graph)
Color the vertices of graph.
Definition GraphColoring.hpp:113