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

Avoid copying the matrix in the ILU preconditioner.

Browse Source 
When the matrix already exists and has the correct sparsity structure
(which is assumed by the calling update()), we can copy just the data,
avoiding the need for allocation and matrix construction.
This commit is contained in:
Atgeirr Flø Rasmussen 2023-02-02 14:10:16 +01:00
parent e6042cdc62
commit 235a8a7f78
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); void reorderBack(const Range& reorderedV, Range& v);
//! \brief The ILU0 decomposition of the matrix. //! \brief The ILU0 decomposition of the matrix.
std::unique_ptr<Matrix> ILU_;
CRS lower_; CRS lower_;
CRS upper_; CRS upper_;
std::vector< block_type > inv_; std::vector< block_type > inv_;

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

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