OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-26 21:36:25 -06:00
Code Issues Packages Projects Releases Wiki Activity

Make ParallelOverlappingILU0 updateable & use it to construct ILU0/n

Browse Source 
This should force the flexible solver approach to recalculate the
decomposition during update.

Closes #2490
This commit is contained in:
Markus Blatt 2020-03-23 09:32:28 +01:00
parent 81b9947e92
commit c14bf078a1
2 changed files with 60 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -21,6 +21,7 @@
#define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
#include <opm/simulators/linalg/GraphColoring.hpp> #include <opm/simulators/linalg/GraphColoring.hpp>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/common/Exceptions.hpp> #include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp> #include <opm/common/ErrorMacros.hpp>
#include <dune/common/version.hh> #include <dune/common/version.hh>
@ -506,6 +507,9 @@ namespace Opm
typedef typename M :: size_type size_type; typedef typename M :: size_type size_type;
lower.clear();
upper.clear();
inv.clear();
lower.resize( A.N() ); lower.resize( A.N() );
upper.resize( A.N() ); upper.resize( A.N() );
inv.resize( A.N() ); inv.resize( A.N() );
@ -599,7 +603,7 @@ namespace Opm
/// used, e.g. Dune::OwnerOverlapCommunication /// used, e.g. Dune::OwnerOverlapCommunication
template<class Matrix, class Domain, class Range, class ParallelInfoT> template<class Matrix, class Domain, class Range, class ParallelInfoT>
class ParallelOverlappingILU0 class ParallelOverlappingILU0
: public Dune::Preconditioner<Domain,Range> : public Dune::PreconditionerWithUpdate<Domain,Range>
{ {
typedef ParallelInfoT ParallelInfo; typedef ParallelInfoT ParallelInfo;
@ -656,6 +660,14 @@ protected:
cols_.push_back( index ); cols_.push_back( index );
} }
void clear()
{
rows_.clear();
values_.clear();
cols_.clear();
nRows_= 0;
}
std::vector< size_type > rows_; std::vector< size_type > rows_;
std::vector< block_type > values_; std::vector< block_type > values_;
std::vector< size_type > cols_; std::vector< size_type > cols_;
@ -691,13 +703,14 @@ public:
upper_(), upper_(),
inv_(), inv_(),
comm_(nullptr), w_(w), comm_(nullptr), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ) relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{ {
interiorSize_ = A.N(); interiorSize_ = A.N();
// BlockMatrix is a Subclass of FieldMatrix that just adds // BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe. // methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), n, milu, redblack, update();
reorder_sphere );
} }
/*! \brief Constructor gets all parameters to operate the prec. /*! \brief Constructor gets all parameters to operate the prec.
@ -721,13 +734,14 @@ public:
upper_(), upper_(),
inv_(), inv_(),
comm_(&comm), w_(w), comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ) relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{ {
interiorSize_ = A.N(); interiorSize_ = A.N();
// BlockMatrix is a Subclass of FieldMatrix that just adds // BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe. // methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), n, milu, redblack, update();
reorder_sphere );
} }
/*! \brief Constructor. /*! \brief Constructor.
@ -772,13 +786,14 @@ public:
upper_(), upper_(),
inv_(), inv_(),
comm_(&comm), w_(w), comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ) relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{ {
interiorSize_ = A.N(); interiorSize_ = A.N();
// BlockMatrix is a Subclass of FieldMatrix that just adds // BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe. // methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), 0, milu, redblack, update();
reorder_sphere );
} }
/*! \brief Constructor. /*! \brief Constructor.
@ -806,12 +821,13 @@ public:
inv_(), inv_(),
comm_(&comm), w_(w), comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ), relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
interiorSize_(interiorSize) interiorSize_(interiorSize),
A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{ {
// BlockMatrix is a Subclass of FieldMatrix that just adds // BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe. // methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), 0, milu, redblack, update( );
reorder_sphere );
} }
/*! /*!
@ -905,15 +921,14 @@ public:
DUNE_UNUSED_PARAMETER(x); DUNE_UNUSED_PARAMETER(x);
} }
protected: virtual void update() override
void init( const Matrix& A, const int iluIteration, MILU_VARIANT milu, bool redBlack, bool reorderSpheres )
{ {
// (For older DUNE versions the communicator might be // (For older DUNE versions the communicator might be
// invalid if redistribution in AMG happened on the coarset level. // invalid if redistribution in AMG happened on the coarset level.
// Therefore we check for nonzero size // Therefore we check for nonzero size
if ( comm_ && comm_->communicator().size()<=0 ) if ( comm_ && comm_->communicator().size()<=0 )
{ {
if ( A.N() > 0 ) if ( A_->N() > 0 )
{ {
OPM_THROW(std::logic_error, "Expected a matrix with zero rows for an invalid communicator."); OPM_THROW(std::logic_error, "Expected a matrix with zero rows for an invalid communicator.");
} }
@ -930,15 +945,15 @@ protected:
std::unique_ptr< Matrix > ILU; std::unique_ptr< Matrix > ILU;
if ( redBlack ) if ( redBlack_ )
{ {
using Graph = Dune::Amg::MatrixGraph<const Matrix>; using Graph = Dune::Amg::MatrixGraph<const Matrix>;
Graph graph(A); Graph graph(*A_);
auto colorsTuple = colorVerticesWelshPowell(graph); auto colorsTuple = colorVerticesWelshPowell(graph);
const auto& colors = std::get<0>(colorsTuple); const auto& colors = std::get<0>(colorsTuple);
const auto& verticesPerColor = std::get<2>(colorsTuple); const auto& verticesPerColor = std::get<2>(colorsTuple);
auto noColors = std::get<1>(colorsTuple); auto noColors = std::get<1>(colorsTuple);
if ( reorderSpheres ) if ( reorderSphere_ )
{ {
ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor, ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor,
graph, 0); graph, 0);
@ -959,27 +974,27 @@ protected:
try try
{ {
if( iluIteration == 0 ) { if( iluIteration_ == 0 ) {
// create ILU-0 decomposition // create ILU-0 decomposition
if ( ordering_.empty() ) if ( ordering_.empty() )
{ {
ILU.reset( new Matrix( A ) ); ILU.reset( new Matrix( *A_ ) );
} }
else else
{ {
ILU.reset( new Matrix(A.N(), A.M(), A.nonzeroes(), Matrix::row_wise)); ILU.reset( new Matrix(A_->N(), A_->M(), 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)
{ {
const auto& row = A[inverseOrdering[iter.index()]]; const auto& row = (*A_)[inverseOrdering[iter.index()]];
for(auto col = row.begin(), cend = row.end(); col != cend; ++col) for(auto col = row.begin(), cend = row.end(); col != cend; ++col)
{ {
iter.insert(ordering_[col.index()]); iter.insert(ordering_[col.index()]);
} }
} }
// Copy values // Copy values
for(auto iter = A.begin(), iend = A.end(); iter != iend; ++iter) for(auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
{ {
auto& newRow = newA[ordering_[iter.index()]]; auto& newRow = newA[ordering_[iter.index()]];
for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col) for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
@ -989,7 +1004,7 @@ protected:
} }
} }
switch ( milu ) switch ( milu_ )
{ {
case MILU_VARIANT::MILU_1: case MILU_VARIANT::MILU_1:
detail::milu0_decomposition ( *ILU); detail::milu0_decomposition ( *ILU);
@ -1007,7 +1022,7 @@ protected:
detail::IsPositiveFunctor() ); detail::IsPositiveFunctor() );
break; break;
default: default:
if (interiorSize_ == A.N()) if (interiorSize_ == A_->N())
bilu0_decomposition( *ILU ); bilu0_decomposition( *ILU );
else else
detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_); detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
@ -1016,7 +1031,7 @@ protected:
} }
else { else {
// create ILU-n decomposition // create ILU-n decomposition
ILU.reset( new Matrix( A.N(), A.M(), Matrix::row_wise) ); ILU.reset( new 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() )
{ {
@ -1029,7 +1044,7 @@ protected:
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)
@ -1053,6 +1068,7 @@ protected:
detail::convertToCRS( *ILU, lower_, upper_, inv_ ); detail::convertToCRS( *ILU, lower_, upper_, inv_ );
} }
protected:
/// \brief Reorder D if needed and return a reference to it. /// \brief Reorder D if needed and return a reference to it.
Range& reorderD(const Range& d) Range& reorderD(const Range& d)
{ {
@ -1123,6 +1139,11 @@ protected:
const field_type w_; const field_type w_;
const bool relaxation_; const bool relaxation_;
size_type interiorSize_; size_type interiorSize_;
const Matrix* A_;
int iluIteration_;
MILU_VARIANT milu_;
bool redBlack_;
bool reorderSphere_;
}; };
} // end namespace Opm } // end namespace Opm

33
opm/simulators/linalg/PreconditionerFactory.hpp
View File

@ -156,26 +156,20 @@ private:
using C = Comm; using C = Comm;
doAddCreator("ILU0", [](const O& op, const P& prm, const C& comm) { doAddCreator("ILU0", [](const O& op, const P& prm, const C& comm) {
const double w = prm.get<double>("relaxation"); const double w = prm.get<double>("relaxation");
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqILU<M, V, V>>>(comm, op.getmat(), w); op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU);
#else
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqILU0<M, V, V>>>(comm, op.getmat(), w);
#endif
}); });
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const C& comm) { doAddCreator("ParOverILU0", [](const O& op, const P& prm, const C& comm) {
const double w = prm.get<double>("relaxation"); const double w = prm.get<double>("relaxation");
// Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner. // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner.
return wrapPreconditioner<Opm::ParallelOverlappingILU0<M, V, V, C>>( return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU); op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU);
}); });
doAddCreator("ILUn", [](const O& op, const P& prm, const C& comm) { doAddCreator("ILUn", [](const O& op, const P& prm, const C& comm) {
const int n = prm.get<int>("ilulevel"); const int n = prm.get<int>("ilulevel");
const double w = prm.get<double>("relaxation"); const double w = prm.get<double>("relaxation");
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqILU<M, V, V>>>(comm, op.getmat(), n, w); op.getmat(), comm, n, w, Opm::MILU_VARIANT::ILU);
#else
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqILUn<M, V, V>>>(comm, op.getmat(), n, w);
#endif
}); });
doAddCreator("Jac", [](const O& op, const P& prm, const C& comm) { doAddCreator("Jac", [](const O& op, const P& prm, const C& comm) {
const int n = prm.get<int>("repeats"); const int n = prm.get<int>("repeats");
@ -229,24 +223,19 @@ private:
using P = boost::property_tree::ptree; using P = boost::property_tree::ptree;
doAddCreator("ILU0", [](const O& op, const P& prm) { doAddCreator("ILU0", [](const O& op, const P& prm) {
const double w = prm.get<double>("relaxation"); const double w = prm.get<double>("relaxation");
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
return wrapPreconditioner<SeqILU<M, V, V>>(op.getmat(), w); op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
#else
return wrapPreconditioner<SeqILU0<M, V, V>>(op.getmat(), w);
#endif
}); });
doAddCreator("ParOverILU0", [](const O& op, const P& prm) { doAddCreator("ParOverILU0", [](const O& op, const P& prm) {
const double w = prm.get<double>("relaxation"); const double w = prm.get<double>("relaxation");
return wrapPreconditioner<Opm::ParallelOverlappingILU0<M, V, V>>(op.getmat(), 0, w, Opm::MILU_VARIANT::ILU); return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
}); });
doAddCreator("ILUn", [](const O& op, const P& prm) { doAddCreator("ILUn", [](const O& op, const P& prm) {
const int n = prm.get<int>("ilulevel"); const int n = prm.get<int>("ilulevel");
const double w = prm.get<double>("relaxation"); const double w = prm.get<double>("relaxation");
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
return wrapPreconditioner<SeqILU<M, V, V>>(op.getmat(), n, w); op.getmat(), n, w, Opm::MILU_VARIANT::ILU);
#else
return wrapPreconditioner<SeqILUn<M, V, V>>(op.getmat(), n, w);
#endif
}); });
doAddCreator("Jac", [](const O& op, const P& prm) { doAddCreator("Jac", [](const O& op, const P& prm) {
const int n = prm.get<int>("repeats"); const int n = prm.get<int>("repeats");