From c14bf078a1159545934f1e41d8e190b5c7f67535 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Mon, 23 Mar 2020 09:32:28 +0100 Subject: [PATCH] Make ParallelOverlappingILU0 updateable & use it to construct ILU0/n This should force the flexible solver approach to recalculate the decomposition during update. Closes #2490 --- .../linalg/ParallelOverlappingILU0.hpp | 77 ++++++++++++------- .../linalg/PreconditionerFactory.hpp | 33 +++----- 2 files changed, 60 insertions(+), 50 deletions(-) diff --git a/opm/simulators/linalg/ParallelOverlappingILU0.hpp b/opm/simulators/linalg/ParallelOverlappingILU0.hpp index 5739fdef0..c5f70e7c1 100644 --- a/opm/simulators/linalg/ParallelOverlappingILU0.hpp +++ b/opm/simulators/linalg/ParallelOverlappingILU0.hpp @@ -21,6 +21,7 @@ #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED #include +#include #include #include #include @@ -506,6 +507,9 @@ namespace Opm typedef typename M :: size_type size_type; + lower.clear(); + upper.clear(); + inv.clear(); lower.resize( A.N() ); upper.resize( A.N() ); inv.resize( A.N() ); @@ -599,7 +603,7 @@ namespace Opm /// used, e.g. Dune::OwnerOverlapCommunication template class ParallelOverlappingILU0 - : public Dune::Preconditioner + : public Dune::PreconditionerWithUpdate { typedef ParallelInfoT ParallelInfo; @@ -656,6 +660,14 @@ protected: cols_.push_back( index ); } + void clear() + { + rows_.clear(); + values_.clear(); + cols_.clear(); + nRows_= 0; + } + std::vector< size_type > rows_; std::vector< block_type > values_; std::vector< size_type > cols_; @@ -691,13 +703,14 @@ public: upper_(), inv_(), comm_(nullptr), w_(w), - relaxation_( std::abs( w - 1.0 ) > 1e-15 ) + 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. - init( reinterpret_cast(A), n, milu, redblack, - reorder_sphere ); + update(); } /*! \brief Constructor gets all parameters to operate the prec. @@ -721,13 +734,14 @@ public: upper_(), inv_(), comm_(&comm), w_(w), - relaxation_( std::abs( w - 1.0 ) > 1e-15 ) + 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. - init( reinterpret_cast(A), n, milu, redblack, - reorder_sphere ); + update(); } /*! \brief Constructor. @@ -772,13 +786,14 @@ public: upper_(), inv_(), comm_(&comm), w_(w), - relaxation_( std::abs( w - 1.0 ) > 1e-15 ) + 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. - init( reinterpret_cast(A), 0, milu, redblack, - reorder_sphere ); + update(); } /*! \brief Constructor. @@ -806,12 +821,13 @@ public: inv_(), comm_(&comm), w_(w), relaxation_( std::abs( w - 1.0 ) > 1e-15 ), - interiorSize_(interiorSize) + 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. - init( reinterpret_cast(A), 0, milu, redblack, - reorder_sphere ); + update( ); } /*! @@ -905,15 +921,14 @@ public: DUNE_UNUSED_PARAMETER(x); } -protected: - void init( const Matrix& A, const int iluIteration, MILU_VARIANT milu, bool redBlack, bool reorderSpheres ) + virtual void update() override { // (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 ) + if ( A_->N() > 0 ) { 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; - if ( redBlack ) + if ( redBlack_ ) { using Graph = Dune::Amg::MatrixGraph; - Graph graph(A); + 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 ( reorderSpheres ) + if ( reorderSphere_ ) { ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor, graph, 0); @@ -959,27 +974,27 @@ protected: try { - if( iluIteration == 0 ) { + if( iluIteration_ == 0 ) { // create ILU-0 decomposition if ( ordering_.empty() ) { - ILU.reset( new Matrix( A ) ); + ILU.reset( new Matrix( *A_ ) ); } 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; // Create sparsity pattern 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) { iter.insert(ordering_[col.index()]); } } // 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()]]; 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: detail::milu0_decomposition ( *ILU); @@ -1007,7 +1022,7 @@ protected: detail::IsPositiveFunctor() ); break; default: - if (interiorSize_ == A.N()) + if (interiorSize_ == A_->N()) bilu0_decomposition( *ILU ); else detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_); @@ -1016,7 +1031,7 @@ protected: } else { // 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 reorderer, inverseReorderer; if ( ordering_.empty() ) { @@ -1029,7 +1044,7 @@ protected: 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) @@ -1053,6 +1068,7 @@ protected: detail::convertToCRS( *ILU, lower_, upper_, inv_ ); } +protected: /// \brief Reorder D if needed and return a reference to it. Range& reorderD(const Range& d) { @@ -1123,6 +1139,11 @@ protected: const field_type w_; const bool relaxation_; size_type interiorSize_; + const Matrix* A_; + int iluIteration_; + MILU_VARIANT milu_; + bool redBlack_; + bool reorderSphere_; }; } // end namespace Opm diff --git a/opm/simulators/linalg/PreconditionerFactory.hpp b/opm/simulators/linalg/PreconditionerFactory.hpp index 0d002e3e3..d7799ebff 100644 --- a/opm/simulators/linalg/PreconditionerFactory.hpp +++ b/opm/simulators/linalg/PreconditionerFactory.hpp @@ -156,26 +156,20 @@ private: using C = Comm; doAddCreator("ILU0", [](const O& op, const P& prm, const C& comm) { const double w = prm.get("relaxation"); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - return wrapBlockPreconditioner>>(comm, op.getmat(), w); -#else - return wrapBlockPreconditioner>>(comm, op.getmat(), w); -#endif + return std::make_shared>( + op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU); }); doAddCreator("ParOverILU0", [](const O& op, const P& prm, const C& comm) { const double w = prm.get("relaxation"); // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner. - return wrapPreconditioner>( + return std::make_shared>( op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU); }); doAddCreator("ILUn", [](const O& op, const P& prm, const C& comm) { const int n = prm.get("ilulevel"); const double w = prm.get("relaxation"); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); -#else - return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); -#endif + return std::make_shared>( + op.getmat(), comm, n, w, Opm::MILU_VARIANT::ILU); }); doAddCreator("Jac", [](const O& op, const P& prm, const C& comm) { const int n = prm.get("repeats"); @@ -229,24 +223,19 @@ private: using P = boost::property_tree::ptree; doAddCreator("ILU0", [](const O& op, const P& prm) { const double w = prm.get("relaxation"); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - return wrapPreconditioner>(op.getmat(), w); -#else - return wrapPreconditioner>(op.getmat(), w); -#endif + return std::make_shared>( + op.getmat(), 0, w, Opm::MILU_VARIANT::ILU); }); doAddCreator("ParOverILU0", [](const O& op, const P& prm) { const double w = prm.get("relaxation"); - return wrapPreconditioner>(op.getmat(), 0, w, Opm::MILU_VARIANT::ILU); + return std::make_shared>( + op.getmat(), 0, w, Opm::MILU_VARIANT::ILU); }); doAddCreator("ILUn", [](const O& op, const P& prm) { const int n = prm.get("ilulevel"); const double w = prm.get("relaxation"); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - return wrapPreconditioner>(op.getmat(), n, w); -#else - return wrapPreconditioner>(op.getmat(), n, w); -#endif + return std::make_shared>( + op.getmat(), n, w, Opm::MILU_VARIANT::ILU); }); doAddCreator("Jac", [](const O& op, const P& prm) { const int n = prm.get("repeats");