OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4427 from atgeirr/do-not-recreate-matrix-ilu0

Browse Source 
Avoid copying the matrix in the ILU preconditioner.
This commit is contained in:
Markus Blatt 2023-02-14 16:11:14 +01:00 committed by GitHub
commit 7449a32d2d
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 31 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
opm/simulators/linalg/ParallelOverlappingILU0.hpp
View File

@ -342,6 +342,7 @@ protected:
void reorderBack(const Range& reorderedV, Range& v);
//! \brief The ILU0 decomposition of the matrix.
std::unique_ptr<Matrix> ILU_;
CRS lower_;
CRS upper_;
std::vector< block_type > inv_;

46
opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp
View File

@ -374,8 +374,6 @@ update()
std::string message;
const int rank = comm_ ? comm_->communicator().rank() : 0;
std::unique_ptr<Matrix> ILU;
if (redBlack_)
{
using Graph = Dune::Amg::MatrixGraph<const Matrix>;
@ -409,13 +407,29 @@ update()
// create ILU-0 decomposition
if (ordering_.empty())
{
ILU = std::make_unique<Matrix>(*A_);
if (ILU_) {
// 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<Matrix>(*A_);
}
}
else
{
ILU = std::make_unique<Matrix>(A_->N(), A_->M(),
A_->nonzeroes(), Matrix::row_wise);
auto& newA = *ILU;
ILU_ = std::make_unique<Matrix>(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)
{
@ -439,35 +453,35 @@ update()
switch (milu_)
{
case MILU_VARIANT::MILU_1:
detail::milu0_decomposition ( *ILU);
detail::milu0_decomposition ( *ILU_);
break;
case MILU_VARIANT::MILU_2:
detail::milu0_decomposition ( *ILU, detail::identityFunctor<typename Matrix::field_type>,
detail::milu0_decomposition ( *ILU_, detail::identityFunctor<typename Matrix::field_type>,
detail::signFunctor<typename Matrix::field_type> );
break;
case MILU_VARIANT::MILU_3:
detail::milu0_decomposition ( *ILU, detail::absFunctor<typename Matrix::field_type>,
detail::milu0_decomposition ( *ILU_, detail::absFunctor<typename Matrix::field_type>,
detail::signFunctor<typename Matrix::field_type> );
break;
case MILU_VARIANT::MILU_4:
detail::milu0_decomposition ( *ILU, detail::identityFunctor<typename Matrix::field_type>,
detail::milu0_decomposition ( *ILU_, detail::identityFunctor<typename Matrix::field_type>,
detail::isPositiveFunctor<typename Matrix::field_type> );
break;
default:
if (interiorSize_ == A_->N())
#if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
bilu0_decomposition( *ILU );
bilu0_decomposition( *ILU_ );
#else
Dune::ILU::blockILU0Decomposition( *ILU );
Dune::ILU::blockILU0Decomposition( *ILU_ );
#endif
else
detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
detail::ghost_last_bilu0_decomposition(*ILU_, interiorSize_);
break;
}
}
else {
// create ILU-n decomposition
ILU = std::make_unique<Matrix>(A_->N(), A_->M(), Matrix::row_wise);
ILU_ = std::make_unique<Matrix>(A_->N(), A_->M(), Matrix::row_wise);
std::unique_ptr<detail::Reorderer> reorderer, inverseReorderer;
if (ordering_.empty())
{
@ -480,7 +494,7 @@ update()
inverseReorderer.reset(new detail::RealReorderer(inverseOrdering));
}
milun_decomposition( *A_, iluIteration_, milu_, *ILU, *reorderer, *inverseReorderer );
milun_decomposition( *A_, iluIteration_, milu_, *ILU_, *reorderer, *inverseReorderer );
}
}
catch (const Dune::MatrixBlockError& error)
@ -501,7 +515,7 @@ update()
}
// store ILU in simple CRS format
detail::convertToCRS(*ILU, lower_, upper_, inv_);
detail::convertToCRS(*ILU_, lower_, upper_, inv_);
}
template<class Matrix, class Domain, class Range, class ParallelInfoT>