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>
43void ghost_last_bilu0_decomposition (M&
A, std::size_t
interiorSize)
49 using block =
typename M::block_type;
59 for (
ij=(*i).begin();
ij.index()<i.index(); ++
ij)
65 (*ij).rightmultiply(*
jj);
72 if (
ik.index()==
jk.index())
81 if (
ik.index()<
jk.index())
89 if (
ij.index()!=i.index())
90 DUNE_THROW(Dune::ISTLError,
"diagonal entry missing");
94 catch (Dune::FMatrixError &
e) {
95 DUNE_THROW(Dune::ISTLError,
"ILU failed to invert matrix block");
101template<
class M,
class CRS,
class InvVector>
102void convertToCRS(
const M&
A, CRS& lower, CRS& upper,
InvVector&
inv)
112 using size_type =
typename M :: size_type;
117 lower.resize(
A.N() );
118 upper.resize(
A.N() );
124 const auto endi =
A.end();
125 for (
auto i =
A.begin(); i !=
endi; ++i) {
126 const size_type
iIndex = i.index();
128 for (
auto j = (*i).begin(); j.index() <
iIndex; ++j) {
136 lower.reserveAdditional(
numLower );
142 for (
auto i=
A.begin(); i!=
endi; ++i, ++
row)
144 const size_type
iIndex = i.index();
147 for (
auto j=(*i).begin(); j.index() <
iIndex; ++j )
149 lower.push_back( (*j), j.index() );
157 const auto rendi =
A.beforeBegin();
162 upper.reserveAdditional(
numUpper );
166 for (
auto i=
A.beforeEnd(); i!=
rendi; --i, ++
row )
168 const size_type
iIndex = i.index();
172 for (
auto j=(*i).beforeEnd(); j.index()>=
iIndex; --j )
174 const size_type
jIndex = j.index();
180 else if ( j.index() >= i.index() )
182 upper.push_back( (*j),
jIndex );
199size_t set_interiorSize(
size_t N,
size_t interiorSize,
const Dune::OwnerOverlapCopyCommunication<int,int>& comm)
208 if (
idx->local().attribute()==1)
210 auto loc =
idx->local().local();
223template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
224Dune::SolverCategory::Category
225ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfoT>::category()
const
227 return std::is_same_v<ParallelInfoT, Dune::Amg::SequentialInformation> ?
228 Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
231template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
241 relaxation_( std::
abs(
w - 1.0 ) > 1
e-15 ),
245 interiorSize_ =
A.N();
251template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
254 const ParallelInfo& comm,
const int n,
const field_type w,
261 relaxation_( std::
abs(
w - 1.0 ) > 1
e-15 ),
265 interiorSize_ =
A.N();
271template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
279template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
289 relaxation_( std::
abs(
w - 1.0 ) > 1
e-15 ),
293 interiorSize_ =
A.N();
299template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
302 const ParallelInfo& comm,
310 relaxation_( std::
abs(
w - 1.0 ) > 1
e-15 ),
321template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
323apply (Domain&
v,
const Range&
d)
326 Range&
md = reorderD(
d);
327 Domain& mv = reorderV(
v);
330 using dblock =
typename Range ::block_type;
331 using vblock =
typename Domain::block_type;
333 const size_type
iEnd = lower_.rows();
337 if (
iEnd != upper_.rows())
339 OPM_THROW(std::logic_error,
"ILU: number of lower and upper rows must be the same");
346 const size_type
rowI = lower_.rows_[ i ];
347 const size_type
rowINext = lower_.rows_[ i+1 ];
351 lower_.values_[
col ].mmv( mv[ lower_.cols_[
col ] ], rhs );
361 const size_type
rowI = upper_.rows_[ i ];
362 const size_type
rowINext = upper_.rows_[ i+1 ];
366 upper_.values_[
col ].mmv( mv[ upper_.cols_[
col ] ], rhs );
370 inv_[ i ].mv( rhs,
vBlock);
373 copyOwnerToAll( mv );
381template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
387 comm_->copyOwnerToAll(
v,
v);
391template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
392void ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfoT>::
399 if (comm_ && comm_->communicator().size() <= 0)
403 OPM_THROW(std::logic_error,
"Expected a matrix with zero rows for an invalid communicator.");
414 const int rank = comm_ ? comm_->communicator().rank() : 0;
418 using Graph = Dune::Amg::MatrixGraph<const Matrix>;
424 if ( reorderSphere_ )
437 std::size_t index = 0;
438 for (
const auto newIndex : ordering_)
446 if (iluIteration_ == 0) {
449 interiorSize_ = detail::set_interiorSize(A_->N(), interiorSize_, *comm_);
454 if (ordering_.empty())
461 for (std::size_t
row = 0;
row < A_->N(); ++
row) {
472 ILU_ = std::make_unique<Matrix>(*A_);
477 ILU_ = std::make_unique<Matrix>(A_->N(), A_->M(),
478 A_->nonzeroes(), Matrix::row_wise);
486 iter.insert(ordering_[
col.index()]);
503 detail::milu0_decomposition ( *ILU_);
506 detail::milu0_decomposition ( *ILU_, detail::identityFunctor<typename Matrix::field_type>,
507 detail::signFunctor<typename Matrix::field_type> );
510 detail::milu0_decomposition ( *ILU_, detail::absFunctor<typename Matrix::field_type>,
511 detail::signFunctor<typename Matrix::field_type> );
514 detail::milu0_decomposition ( *ILU_, detail::identityFunctor<typename Matrix::field_type>,
515 detail::isPositiveFunctor<typename Matrix::field_type> );
518 if (interiorSize_ == A_->N())
519 Dune::ILU::blockILU0Decomposition( *ILU_ );
521 detail::ghost_last_bilu0_decomposition(*ILU_, interiorSize_);
527 ILU_ = std::make_unique<Matrix>(A_->N(), A_->M(), Matrix::row_wise);
529 if (ordering_.empty())
531 reorderer.reset(
new detail::NoReorderer());
536 reorderer.reset(
new detail::RealReorderer(ordering_));
543 catch (
const Dune::MatrixBlockError& error)
545 message = error.what();
546 std::cerr <<
"Exception occurred on process " << rank <<
" during " <<
547 "setup of ILU0 preconditioner with message: "
548 << message<<std::endl;
557 throw Dune::MatrixBlockError();
561 detail::convertToCRS(*ILU_, lower_, upper_, inv_);
564template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
568 if (ordering_.empty())
574 return const_cast<Range&
>(
d);
578 reorderedD_.resize(
d.size());
580 for (
const auto index : ordering_)
582 reorderedD_[index] =
d[i++];
588template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
592 if (ordering_.empty())
598 reorderedV_.resize(
v.size());
600 for (
const auto index : ordering_)
602 reorderedV_[index] =
v[i++];
608template<
class Matrix,
class Domain,
class Range,
class ParallelInfoT>
612 if (!ordering_.empty())
615 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:233
Domain & reorderV(Domain &v)
Reorder V if needed and return a reference to it.
Definition ParallelOverlappingILU0_impl.hpp:590
Range & reorderD(const Range &d)
Reorder D if needed and return a reference to it.
Definition ParallelOverlappingILU0_impl.hpp:566
void apply(Domain &v, const Range &d) override
Apply the preconditoner.
Definition ParallelOverlappingILU0_impl.hpp:323
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