From 6d21214fa709aff9ce9ef700fe508a7dff068f5b Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 1 Feb 2018 08:50:16 +0100 Subject: [PATCH] Switch to ParallelOverlappingILU0 for CPRPreconditioner. This seems to have been forgotten previously. Now the code int CPRPreconditioner.hpp uses ParallelOverlappingILU0 instead of SeqILU[0n]/BlockPreconditioner which makes the code more slim. --- CMakeLists_files.cmake | 1 - opm/autodiff/AdditionalObjectDeleter.hpp | 58 ----- opm/autodiff/CPRPreconditioner.hpp | 206 ++++-------------- opm/autodiff/ISTLSolver.hpp | 3 +- .../NewtonIterationBlackoilInterleaved.cpp | 1 - opm/autodiff/ParallelOverlappingILU0.hpp | 20 ++ 6 files changed, 63 insertions(+), 226 deletions(-) delete mode 100644 opm/autodiff/AdditionalObjectDeleter.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 5f53b3021..3714c2195 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -252,7 +252,6 @@ list (APPEND PROGRAM_SOURCE_FILES # originally generated with the command: # find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort list (APPEND PUBLIC_HEADER_FILES - opm/autodiff/AdditionalObjectDeleter.hpp opm/autodiff/AutoDiffBlock.hpp opm/autodiff/AutoDiffHelpers.hpp opm/autodiff/AutoDiffMatrix.hpp diff --git a/opm/autodiff/AdditionalObjectDeleter.hpp b/opm/autodiff/AdditionalObjectDeleter.hpp deleted file mode 100644 index 985bf741f..000000000 --- a/opm/autodiff/AdditionalObjectDeleter.hpp +++ /dev/null @@ -1,58 +0,0 @@ - /* - Copyright 2015 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 . -*/ -#ifndef OPM_ADDITIONALOBJECTDELETER_HEADER_INCLUDED -#define OPM_ADDITIONALOBJECTDELETER_HEADER_INCLUDED -namespace Opm -{ -//! \brief A custom deleter that will delete an additional object, too. -//! -//! In dune-istl most parallel preconditioners hold a reference to -//! a sequential preconditioner. -//! In CPRPreconditioner and NewtonIterationBlackoilInterleaved we use unique_ptr -//! for the memory management. -//! Ergo we need to construct the sequential preconditioner with new and -//! make sure that it gets deleted together with the enclosing parallel -//! preconditioner. Therefore this deleter stores a pointer to it and deletes -//! it during destruction. -//! \tparam The type of the additional object to be deleted. -template -class AdditionalObjectDeleter -{ -public: - //! \brief empty constructor. - AdditionalObjectDeleter() - : additional_object_() - {} - //! \brief Constructor taking the object that needs to deleted. - AdditionalObjectDeleter(T& additional_object) - : additional_object_(&additional_object){} - //! \brief Delete an object and the additional one. - template - void operator()(T1* pt) - { - delete pt; - delete additional_object_; - } -private: - T* additional_object_; -}; - -} -#endif diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index 9a0862d5b..de1f92ded 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -47,7 +47,6 @@ #include #include -#include #include #include namespace Opm @@ -67,47 +66,13 @@ template struct CPRSelector { /// \brief The information about the parallelization and communication - typedef Dune::Amg::SequentialInformation ParallelInformation; - /// \brief The operator type; - typedef Dune::MatrixAdapter Operator; - /// \brief The type of the preconditioner used for the elliptic part. - typedef Dune::SeqILU0 EllipticPreconditioner; - /// \brief The type of the unique pointer to the preconditioner of the elliptic part. - typedef std::unique_ptr EllipticPreconditionerPointer; - - /// \brief type of AMG used to precondition the elliptic system. - typedef EllipticPreconditioner Smoother; - typedef Dune::Amg::AMG AMG; - - /// \brief creates an Operator from the matrix - /// \param M The matrix to use. - /// \param p The parallel information to use. - static Operator* makeOperator(const M& m, const P&) - { - return new Operator(m); - } - -}; - -#if HAVE_MPI -/// \copydoc CPRSelector -template -struct CPRSelector > -{ - /// \brief The information about the parallelization and communication - typedef Dune::OwnerOverlapCopyCommunication ParallelInformation; + using ParallelInformation = P; /// \brief The operator type; typedef Dune::OverlappingSchwarzOperator Operator; - /// \brief The type of the preconditioner used for the elliptic part. - //typedef Dune::BlockPreconditioner > - //EllipticPreconditioner; - //typedef ParallelRestrictedOverlappingSchwarz > - //EllipticPreconditioner; typedef ParallelOverlappingILU0 EllipticPreconditioner; /// \brief The type of the unique pointer to the preconditioner of the elliptic part. - typedef std::unique_ptr > > + typedef std::unique_ptr EllipticPreconditionerPointer; typedef EllipticPreconditioner Smoother; @@ -122,155 +87,68 @@ struct CPRSelector > } }; -//! \brief Creates the deleter needed for the parallel ILU preconditioners. -//! \tparam ILU The type of the underlying sequential ILU preconditioner. -//! \tparam I1 The global index type. -//! \tparam I2 The local index type. -//! \param ilu A reference to the wrapped preconditioner -//! \param p The parallel information for template parameter deduction. -template -AdditionalObjectDeleter -createParallelDeleter(ILU& ilu, const Dune::OwnerOverlapCopyCommunication& p) - { - (void) p; - return AdditionalObjectDeleter(ilu); - } -#endif +template +struct CPRSelector +{ + /// \brief The information about the parallelization and communication + typedef Dune::Amg::SequentialInformation ParallelInformation; + /// \brief The operator type; + typedef Dune::MatrixAdapter Operator; + /// \brief The type of the preconditioner used for the elliptic part. + typedef ParallelOverlappingILU0 + EllipticPreconditioner; + /// \brief The type of the unique pointer to the preconditioner of the elliptic part. + typedef std::unique_ptr EllipticPreconditionerPointer; + /// \brief type of AMG used to precondition the elliptic system. + typedef EllipticPreconditioner Smoother; + typedef Dune::Amg::AMG AMG; + + /// \brief creates an Operator from the matrix + /// \param M The matrix to use. + /// \param p The parallel information to use. + static Operator* makeOperator(const M& m, const ParallelInformation&) + { + return new Operator(m); + } + +}; //! \brief Creates and initializes a unique pointer to an sequential ILU0 preconditioner. //! \param A The matrix of the linear system to solve. //! \param relax The relaxation factor to use. -template -std::shared_ptr > -createILU0Ptr(const M& A, double relax, const Dune::Amg::SequentialInformation&) +template +std::shared_ptr > +createILU0Ptr(const M& A, const C& comm, double relax) { - return std::shared_ptr >(new Dune::SeqILU0( A, relax) ); + return std::make_shared >(A, comm, relax); } //! \brief Creates and initializes a shared pointer to an ILUn preconditioner. //! \param A The matrix of the linear system to solve. //! \param ilu_n The n parameter for the extension of the nonzero pattern. //! \param relax The relaxation factor to use. -template -std::shared_ptr > -createILUnPtr(const M& A, int ilu_n, double relax, const Dune::Amg::SequentialInformation&) +template +std::shared_ptr > +createILUnPtr(const M& A, const C& comm, int ilu_n, double relax) { - return std::shared_ptr >(new Dune::SeqILUn( A, ilu_n, relax) ); + return std::make_shared >( A, comm, ilu_n, relax); } - -#if HAVE_MPI -template -struct SelectParallelILUSharedPtr -{ - typedef std::shared_ptr< - Dune::BlockPreconditioner< - typename ILU::range_type, - typename ILU::domain_type, - Dune::OwnerOverlapCopyCommunication, - ILU - > - > type; -}; - -//! \brief Creates and initializes a shared pointer to an ILUn preconditioner. -//! \param A The matrix of the linear system to solve. -//! \param relax The relaxation factor to use. -/// \param comm The object describing the parallelization information and communication. -template -typename SelectParallelILUSharedPtr, I1, I2>::type -createILU0Ptr(const M& A, double relax, - const Dune::OwnerOverlapCopyCommunication& comm) -{ - typedef Dune::BlockPreconditioner< - X, - X, - Dune::OwnerOverlapCopyCommunication, - Dune::SeqILU0 - > PointerType; - Dune::SeqILU0* ilu = nullptr; - int ilu_setup_successful = 1; - std::string message; - try { - ilu = new Dune::SeqILU0(A, relax); - } - catch ( Dune::MatrixBlockError error ) - { - message = error.what(); - std::cerr<<"Exception occured on process " << - comm.communicator().rank() << " during " << - "setup of ILU0 preconditioner with message: " << - message<, I1, I2> - ::type ( new PointerType(*ilu, comm), createParallelDeleter(*ilu, comm)); -} - -//! \brief Creates and initializes a shared pointer to an ILUn preconditioner. -//! \param A The matrix of the linear system to solve. -//! \param ilu_n The n parameter for the extension of the nonzero pattern. -//! \param relax The relaxation factor to use. -/// \param comm The object describing the parallelization information and communication. -template -typename SelectParallelILUSharedPtr, I1, I2>::type -createILUnPtr(const M& A, int ilu_n, double relax, - const Dune::OwnerOverlapCopyCommunication& comm) -{ - typedef Dune::BlockPreconditioner< - X, - X, - Dune::OwnerOverlapCopyCommunication, - Dune::SeqILUn - > PointerType; - Dune::SeqILUn* ilu = new Dune::SeqILUn( A, ilu_n, relax); - - return typename SelectParallelILUSharedPtr, I1, I2>::type - (new PointerType(*ilu, comm),createParallelDeleter(*ilu, comm)); -} -#endif - -/// \brief Creates the elliptic preconditioner (ILU0) -/// \param Ae The matrix of the elliptic system. -/// \param relax The relaxation parameter for ILU0 -template -std::unique_ptr > -createEllipticPreconditionerPointer(const M& Ae, double relax, - const Dune::Amg::SequentialInformation&) -{ - return std::unique_ptr >(new Dune::SeqILU0(Ae, relax)); -} - -#if HAVE_MPI /// \brief Creates the elliptic preconditioner (ILU0) /// \param Ae The matrix of the elliptic system. /// \param relax The relaxation parameter for ILU0. /// \param comm The object describing the parallelization information and communication. -template -typename CPRSelector > -::EllipticPreconditionerPointer +template +typename CPRSelector::EllipticPreconditionerPointer createEllipticPreconditionerPointer(const M& Ae, double relax, - const Dune::OwnerOverlapCopyCommunication& comm) + const P& comm) { - typedef typename CPRSelector > + typedef typename CPRSelector ::EllipticPreconditioner ParallelPreconditioner; - //Dune::SeqILU0* ilu=new Dune::SeqILU0(Ae, relax); - typedef typename CPRSelector > + typedef typename CPRSelector ::EllipticPreconditionerPointer EllipticPreconditionerPointer; return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax)); } -#endif // #if HAVE_MPI - template < class C, class Op, class P, class AMG > inline void createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr ) @@ -456,10 +334,10 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std: // create the preconditioner for the whole system. if( param_.cpr_ilu_n_ == 0 ) { - pre_ = ISTLUtility::createILU0Ptr( A_, param_.cpr_relax_, comm_ ); + pre_ = ISTLUtility::createILU0Ptr( A_, comm, param_.cpr_relax_ ); } else { - pre_ = ISTLUtility::createILUnPtr( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm_ ); + pre_ = ISTLUtility::createILUnPtr( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_); } } diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp index 100ab073b..cc8089e05 100644 --- a/opm/autodiff/ISTLSolver.hpp +++ b/opm/autodiff/ISTLSolver.hpp @@ -20,7 +20,6 @@ #ifndef OPM_ISTLSOLVER_HEADER_INCLUDED #define OPM_ISTLSOLVER_HEADER_INCLUDED -#include #include #include #include @@ -418,7 +417,7 @@ namespace Opm opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) ); } - const double relax = 1.0; + const double relax = parameters_.ilu_relaxation_; if ( ! parameters_.amg_blackoil_system_ ) { typedef typename CPRSelectorType::AMG AMG; diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 7db49042c..afaa217f2 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -27,7 +27,6 @@ #define FLOW_SUPPORT_AMG !defined(HAVE_UMFPACK) #include -#include #include #include #include diff --git a/opm/autodiff/ParallelOverlappingILU0.hpp b/opm/autodiff/ParallelOverlappingILU0.hpp index bfa7fb2c1..f72affd66 100644 --- a/opm/autodiff/ParallelOverlappingILU0.hpp +++ b/opm/autodiff/ParallelOverlappingILU0.hpp @@ -265,6 +265,26 @@ public: init( reinterpret_cast(A), n ); } + /*! \brief Constructor gets all parameters to operate the prec. + \param A The matrix to operate on. + \param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication + \param n ILU fill in level (for testing). This does not work in parallel. + \param w The relaxation factor. + */ + template + ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, + const ParallelInfo& comm, const int n, const field_type w ) + : lower_(), + upper_(), + inv_(), + comm_(&comm), w_(w), + relaxation_( std::abs( w - 1.0 ) > 1e-15 ) + { + // BlockMatrix is a Subclass of FieldMatrix that just adds + // methods. Therefore this cast should be safe. + init( reinterpret_cast(A), n ); + } + /*! \brief Constructor. Constructor gets all parameters to operate the prec.