/*
Copyright 2015, 2022 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 Statoil AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see .
*/
#include
#include
#include
#include
#include
#include
#include
#include
namespace Opm
{
namespace detail
{
//! Compute Blocked ILU0 decomposition, when we know junk ghost rows are located at the end of A
template
void ghost_last_bilu0_decomposition (M& A, size_t interiorSize)
{
// iterator types
using rowiterator = typename M::RowIterator;
using coliterator = typename M::ColIterator;
using block = typename M::block_type;
// implement left looking variant with stored inverse
for (rowiterator i = A.begin(); i.index() < interiorSize; ++i)
{
// coliterator is diagonal after the following loop
coliterator endij=(*i).end(); // end of row i
coliterator ij;
// eliminate entries left of diagonal; store L factor
for (ij=(*i).begin(); ij.index()
void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv)
{
OPM_TIMEBLOCK(convertToCRS);
// No need to do anything for 0 rows. Return to prevent indexing a
// a zero sized array.
if ( A.N() == 0 )
{
return;
}
using size_type = typename M :: size_type;
lower.clear();
upper.clear();
inv.clear();
lower.resize( A.N() );
upper.resize( A.N() );
inv.resize( A.N() );
// Count the lower and upper matrix entries.
size_type numLower = 0;
size_type numUpper = 0;
const auto endi = A.end();
for (auto i = A.begin(); i != endi; ++i) {
const size_type iIndex = i.index();
size_type numLowerRow = 0;
for (auto j = (*i).begin(); j.index() < iIndex; ++j) {
++numLowerRow;
}
numLower += numLowerRow;
numUpper += (*i).size() - numLowerRow - 1;
}
assert(numLower + numUpper + A.N() == A.nonzeroes());
lower.reserveAdditional( numLower );
// implement left looking variant with stored inverse
size_type row = 0;
size_type colcount = 0;
lower.rows_[ 0 ] = colcount;
for (auto i=A.begin(); i!=endi; ++i, ++row)
{
const size_type iIndex = i.index();
// eliminate entries left of diagonal; store L factor
for (auto j=(*i).begin(); j.index() < iIndex; ++j )
{
lower.push_back( (*j), j.index() );
++colcount;
}
lower.rows_[ iIndex+1 ] = colcount;
}
assert(colcount == numLower);
const auto rendi = A.beforeBegin();
row = 0;
colcount = 0;
upper.rows_[ 0 ] = colcount ;
upper.reserveAdditional( numUpper );
// NOTE: upper and inv store entries in reverse order, reverse here
// relative to ILU
for (auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
{
const size_type iIndex = i.index();
// store in reverse row order
// eliminate entries left of diagonal; store L factor
for (auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
{
const size_type jIndex = j.index();
if( j.index() == iIndex )
{
inv[ row ] = (*j);
break;
}
else if ( j.index() >= i.index() )
{
upper.push_back( (*j), jIndex );
++colcount ;
}
}
upper.rows_[ row+1 ] = colcount;
}
assert(colcount == numUpper);
}
} // end namespace detail
template
Dune::SolverCategory::Category
ParallelOverlappingILU0::category() const
{
return std::is_same_v ?
Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
}
template
ParallelOverlappingILU0::
ParallelOverlappingILU0(const Matrix& A,
const int n, const field_type w,
MILU_VARIANT milu, bool redblack,
bool reorder_sphere)
: lower_(),
upper_(),
inv_(),
comm_(nullptr), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
A_(&reinterpret_cast(A)), iluIteration_(n),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{
interiorSize_ = A.N();
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
update();
}
template
ParallelOverlappingILU0::
ParallelOverlappingILU0(const Matrix& A,
const ParallelInfo& comm, const int n, const field_type w,
MILU_VARIANT milu, bool redblack,
bool reorder_sphere)
: lower_(),
upper_(),
inv_(),
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
A_(&reinterpret_cast(A)), iluIteration_(n),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{
interiorSize_ = A.N();
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
update();
}
template
ParallelOverlappingILU0::
ParallelOverlappingILU0(const Matrix& A,
const field_type w, MILU_VARIANT milu, bool redblack,
bool reorder_sphere)
: ParallelOverlappingILU0( A, 0, w, milu, redblack, reorder_sphere )
{}
template
ParallelOverlappingILU0::
ParallelOverlappingILU0(const Matrix& A,
const ParallelInfo& comm, const field_type w,
MILU_VARIANT milu, bool redblack,
bool reorder_sphere)
: lower_(),
upper_(),
inv_(),
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
A_(&reinterpret_cast(A)), iluIteration_(0),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{
interiorSize_ = A.N();
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
update();
}
template
ParallelOverlappingILU0::
ParallelOverlappingILU0(const Matrix& A,
const ParallelInfo& comm,
const field_type w, MILU_VARIANT milu,
size_type interiorSize, bool redblack,
bool reorder_sphere)
: lower_(),
upper_(),
inv_(),
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
interiorSize_(interiorSize),
A_(&reinterpret_cast(A)), iluIteration_(0),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
update( );
}
template
void ParallelOverlappingILU0::
apply (Domain& v, const Range& d)
{
OPM_TIMEBLOCK(apply);
Range& md = reorderD(d);
Domain& mv = reorderV(v);
// iterator types
using dblock = typename Range ::block_type;
using vblock = typename Domain::block_type;
const size_type iEnd = lower_.rows();
const size_type lastRow = iEnd - 1;
size_type upperLoopStart = iEnd - interiorSize_;
size_type lowerLoopEnd = interiorSize_;
if (iEnd != upper_.rows())
{
OPM_THROW(std::logic_error,"ILU: number of lower and upper rows must be the same");
}
// lower triangular solve
for (size_type i = 0; i < lowerLoopEnd; ++i)
{
dblock rhs( md[ i ] );
const size_type rowI = lower_.rows_[ i ];
const size_type rowINext = lower_.rows_[ i+1 ];
for (size_type col = rowI; col < rowINext; ++col)
{
lower_.values_[ col ].mmv( mv[ lower_.cols_[ col ] ], rhs );
}
mv[ i ] = rhs; // Lii = I
}
for (size_type i = upperLoopStart; i < iEnd; ++i)
{
vblock& vBlock = mv[ lastRow - i ];
vblock rhs ( vBlock );
const size_type rowI = upper_.rows_[ i ];
const size_type rowINext = upper_.rows_[ i+1 ];
for (size_type col = rowI; col < rowINext; ++col)
{
upper_.values_[ col ].mmv( mv[ upper_.cols_[ col ] ], rhs );
}
// apply inverse and store result
inv_[ i ].mv( rhs, vBlock);
}
copyOwnerToAll( mv );
if( relaxation_ ) {
mv *= w_;
}
reorderBack(mv, v);
}
template
template
void ParallelOverlappingILU0::
copyOwnerToAll(V& v) const
{
if( comm_ ) {
comm_->copyOwnerToAll(v, v);
}
}
template
void ParallelOverlappingILU0::
update()
{
OPM_TIMEBLOCK(update);
// (For older DUNE versions the communicator might be
// invalid if redistribution in AMG happened on the coarset level.
// Therefore we check for nonzero size
if (comm_ && comm_->communicator().size() <= 0)
{
if (A_->N() > 0)
{
OPM_THROW(std::logic_error, "Expected a matrix with zero rows for an invalid communicator.");
}
else
{
// simply set the communicator to null
comm_ = nullptr;
}
}
int ilu_setup_successful = 1;
std::string message;
const int rank = comm_ ? comm_->communicator().rank() : 0;
if (redBlack_)
{
using Graph = Dune::Amg::MatrixGraph;
Graph graph(*A_);
auto colorsTuple = colorVerticesWelshPowell(graph);
const auto& colors = std::get<0>(colorsTuple);
const auto& verticesPerColor = std::get<2>(colorsTuple);
auto noColors = std::get<1>(colorsTuple);
if ( reorderSphere_ )
{
ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor,
graph, 0);
}
else
{
ordering_ = reorderVerticesPreserving(colors, noColors, verticesPerColor,
graph);
}
}
std::vector inverseOrdering(ordering_.size());
std::size_t index = 0;
for (const auto newIndex : ordering_)
{
inverseOrdering[newIndex] = index++;
}
try
{
OPM_TIMEBLOCK(iluDecomposition);
if (iluIteration_ == 0) {
// create ILU-0 decomposition
if (ordering_.empty())
{
if (ILU_) {
OPM_TIMEBLOCK(iluDecompositionMakeMatrix);
// The ILU_ matrix is already a copy with the same
// sparse structure as A_, but the values of A_ may
// have changed, so we must copy all elements.
for (size_t row = 0; row < A_->N(); ++row) {
const auto& Arow = (*A_)[row];
auto& ILUrow = (*ILU_)[row];
auto Ait = Arow.begin();
auto Iit = ILUrow.begin();
for (; Ait != Arow.end(); ++Ait, ++Iit) {
*Iit = *Ait;
}
}
} else {
// First call, must duplicate matrix.
ILU_ = std::make_unique(*A_);
}
}
else
{
ILU_ = std::make_unique(A_->N(), A_->M(),
A_->nonzeroes(), Matrix::row_wise);
auto& newA = *ILU_;
// Create sparsity pattern
for (auto iter = newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
{
const auto& row = (*A_)[inverseOrdering[iter.index()]];
for (auto col = row.begin(), cend = row.end(); col != cend; ++col)
{
iter.insert(ordering_[col.index()]);
}
}
// Copy values
for (auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
{
auto& newRow = newA[ordering_[iter.index()]];
for (auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
{
newRow[ordering_[col.index()]] = *col;
}
}
}
switch (milu_)
{
case MILU_VARIANT::MILU_1:
detail::milu0_decomposition ( *ILU_);
break;
case MILU_VARIANT::MILU_2:
detail::milu0_decomposition ( *ILU_, detail::identityFunctor,
detail::signFunctor );
break;
case MILU_VARIANT::MILU_3:
detail::milu0_decomposition ( *ILU_, detail::absFunctor,
detail::signFunctor );
break;
case MILU_VARIANT::MILU_4:
detail::milu0_decomposition ( *ILU_, detail::identityFunctor,
detail::isPositiveFunctor );
break;
default:
if (interiorSize_ == A_->N())
#if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
bilu0_decomposition( *ILU_ );
#else
Dune::ILU::blockILU0Decomposition( *ILU_ );
#endif
else
detail::ghost_last_bilu0_decomposition(*ILU_, interiorSize_);
break;
}
}
else {
// create ILU-n decomposition
ILU_ = std::make_unique(A_->N(), A_->M(), Matrix::row_wise);
std::unique_ptr reorderer, inverseReorderer;
if (ordering_.empty())
{
reorderer.reset(new detail::NoReorderer());
inverseReorderer.reset(new detail::NoReorderer());
}
else
{
reorderer.reset(new detail::RealReorderer(ordering_));
inverseReorderer.reset(new detail::RealReorderer(inverseOrdering));
}
milun_decomposition( *A_, iluIteration_, milu_, *ILU_, *reorderer, *inverseReorderer );
}
}
catch (const Dune::MatrixBlockError& error)
{
message = error.what();
std::cerr << "Exception occurred on process " << rank << " during " <<
"setup of ILU0 preconditioner with message: "
<< message<communicator().min(ilu_setup_successful) == 0;
const bool local_failure = ilu_setup_successful == 0;
if (local_failure || parallel_failure)
{
throw Dune::MatrixBlockError();
}
// store ILU in simple CRS format
detail::convertToCRS(*ILU_, lower_, upper_, inv_);
}
template
Range& ParallelOverlappingILU0::
reorderD(const Range& d)
{
if (ordering_.empty())
{
// As d is non-const in the apply method of the
// solver casting away constness in this particular
// setting is not undefined. It is ugly though but due
// to the preconditioner interface of dune-istl.
return const_cast(d);
}
else
{
reorderedD_.resize(d.size());
std::size_t i = 0;
for (const auto index : ordering_)
{
reorderedD_[index] = d[i++];
}
return reorderedD_;
}
}
template
Domain& ParallelOverlappingILU0::
reorderV(Domain& v)
{
if (ordering_.empty())
{
return v;
}
else
{
reorderedV_.resize(v.size());
std::size_t i = 0;
for (const auto index : ordering_)
{
reorderedV_[index] = v[i++];
}
return reorderedV_;
}
}
template
void ParallelOverlappingILU0::
reorderBack(const Range& reorderedV, Range& v)
{
if (!ordering_.empty())
{
std::size_t i = 0;
for (const auto index : ordering_)
{
v[i++] = reorderedV[index];
}
}
}
} // end namespace Opm