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

@@ -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

opm/autodiff/AdditionalObjectDeleter.hpp

@@ -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 <http://www.gnu.org/licenses/>.
-*/
-#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 T>
-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<class T1>
-    void operator()(T1* pt)
-    {
-        delete pt;
-        delete additional_object_;
-    }
-private:
-    T* additional_object_;
-};
-
-}
-#endif

opm/autodiff/CPRPreconditioner.hpp

@@ -47,7 +47,6 @@
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/Exceptions.hpp>
 
-#include <opm/autodiff/AdditionalObjectDeleter.hpp>
 #include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
 #include <opm/autodiff/ParallelOverlappingILU0.hpp>
 namespace Opm
@@ -67,47 +66,13 @@ template<class M, class X, class Y, class P>
 struct CPRSelector
 {
     /// \brief The information about the parallelization and communication
-    typedef Dune::Amg::SequentialInformation ParallelInformation;
+    using ParallelInformation = P;
-    /// \brief The operator type;
-    typedef Dune::MatrixAdapter<M, X, Y> Operator;
-    /// \brief The type of the preconditioner used for the elliptic part.
-    typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
-    /// \brief The type of the unique pointer to the preconditioner of the elliptic part.
-    typedef std::unique_ptr<EllipticPreconditioner> EllipticPreconditionerPointer;
-
-    /// \brief type of AMG used to precondition the elliptic system.
-    typedef EllipticPreconditioner Smoother;
-    typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> 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<M,X,X,Y,P>
-template<class M, class X, class Y, class I1, class I2>
-struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
-{
-    /// \brief The information about the parallelization and communication
-    typedef Dune::OwnerOverlapCopyCommunication<I1,I2> ParallelInformation;
     /// \brief The operator type;
     typedef Dune::OverlappingSchwarzOperator<M,X,X,ParallelInformation> Operator;
-    /// \brief The type of the preconditioner used for the elliptic part.
-    //typedef Dune::BlockPreconditioner<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
-    //EllipticPreconditioner;
-    //typedef ParallelRestrictedOverlappingSchwarz<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
-    //EllipticPreconditioner;
     typedef ParallelOverlappingILU0<M,X, X, ParallelInformation>
     EllipticPreconditioner;
     /// \brief The type of the unique pointer to the preconditioner of the elliptic part.
-    typedef std::unique_ptr<EllipticPreconditioner,
+    typedef std::unique_ptr<EllipticPreconditioner>
-                            AdditionalObjectDeleter<Dune::SeqILU0<M,X,X> > >
     EllipticPreconditionerPointer;
 
     typedef EllipticPreconditioner Smoother;
@@ -122,155 +87,68 @@ struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
     }
 };
 
-//! \brief Creates the deleter needed for the parallel ILU preconditioners.
+template<class M, class X, class Y>
-//! \tparam ILU The type of the underlying sequential ILU preconditioner.
+struct CPRSelector<M,X,Y,Dune::Amg::SequentialInformation>
-//! \tparam I1 The global index type.
+{
-//! \tparam I2 The local index type.
+    /// \brief The information about the parallelization and communication
-//! \param  ilu A reference to the wrapped preconditioner
+    typedef Dune::Amg::SequentialInformation ParallelInformation;
-//! \param  p The parallel information for template parameter deduction.
+    /// \brief The operator type;
-template<class ILU, class I1, class I2>
+    typedef Dune::MatrixAdapter<M, X, Y> Operator;
-AdditionalObjectDeleter<ILU>
+    /// \brief The type of the preconditioner used for the elliptic part.
-createParallelDeleter(ILU& ilu, const Dune::OwnerOverlapCopyCommunication<I1,I2>& p)
+    typedef ParallelOverlappingILU0<M,X, X, ParallelInformation>
-    {
+    EllipticPreconditioner;
-        (void) p;
+    /// \brief The type of the unique pointer to the preconditioner of the elliptic part.
-        return AdditionalObjectDeleter<ILU>(ilu);
+    typedef std::unique_ptr<EllipticPreconditioner> EllipticPreconditionerPointer;
-    }
-#endif
 
+    /// \brief type of AMG used to precondition the elliptic system.
+    typedef EllipticPreconditioner Smoother;
+    typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> 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<class M, class X>
+template<class M, class X, class C>
-std::shared_ptr<Dune::SeqILU0<M,X,X> >
+std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
-createILU0Ptr(const M& A, double relax, const Dune::Amg::SequentialInformation&)
+createILU0Ptr(const M& A, const C& comm, double relax)
 {
-    return std::shared_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>( A, relax) );
+    return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >(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<class M, class X>
+template<class M, class X, class C>
-std::shared_ptr<Dune::SeqILUn<M,X,X> >
+std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
-createILUnPtr(const M& A, int ilu_n, double relax, const Dune::Amg::SequentialInformation&)
+createILUnPtr(const M& A, const C& comm, int ilu_n, double relax)
 {
-    return std::shared_ptr<Dune::SeqILUn<M,X,X> >(new Dune::SeqILUn<M,X,X>( A, ilu_n, relax) );
+    return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >( A, comm, ilu_n, relax);
 }
-
-#if HAVE_MPI
-template<class ILU, class I1, class I2>
-struct SelectParallelILUSharedPtr
-{
-    typedef std::shared_ptr<
-        Dune::BlockPreconditioner<
-            typename ILU::range_type,
-            typename ILU::domain_type,
-            Dune::OwnerOverlapCopyCommunication<I1,I2>,
-            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<class M, class X, class I1, class I2>
-typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>::type
-createILU0Ptr(const M& A, double relax,
-              const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
-{
-    typedef Dune::BlockPreconditioner<
-        X,
-        X,
-        Dune::OwnerOverlapCopyCommunication<I1,I2>,
-        Dune::SeqILU0<M,X,X>
-        > PointerType;
-    Dune::SeqILU0<M,X,X>* ilu = nullptr;
-    int ilu_setup_successful = 1;
-    std::string message;
-    try {
-        ilu = new Dune::SeqILU0<M,X,X>(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<<std::endl;
-                ilu_setup_successful = 0;
-    }
-    // Check whether there was a problem on some process
-    if ( comm.communicator().min(ilu_setup_successful) == 0 )
-    {
-        if ( ilu ) // not null if constructor succeeded
-        {
-            // prevent memory leak
-            delete ilu;
-        }
-        throw Dune::MatrixBlockError();
-    }
-    return typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, 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<class M, class X, class I1, class I2>
-typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
-createILUnPtr(const M& A, int ilu_n, double relax,
-              const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
-{
-    typedef Dune::BlockPreconditioner<
-        X,
-        X,
-        Dune::OwnerOverlapCopyCommunication<I1,I2>,
-        Dune::SeqILUn<M,X,X>
-        > PointerType;
-    Dune::SeqILUn<M,X,X>* ilu = new Dune::SeqILUn<M,X,X>( A, ilu_n, relax);
-
-    return typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, 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<class M, class X=typename M::range_type>
-std::unique_ptr<Dune::SeqILU0<M,X,X> >
-createEllipticPreconditionerPointer(const M& Ae, double relax,
-                                    const Dune::Amg::SequentialInformation&)
-{
-    return std::unique_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>(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<class M, class X=typename M::range_type, class I1, class I2>
+template<class M, class X=typename M::range_type, class P>
-typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
+typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
-::EllipticPreconditionerPointer
 createEllipticPreconditionerPointer(const M& Ae, double relax,
-                                    const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
+                                    const P& comm)
 {
-    typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
+    typedef typename CPRSelector<M,X,X,P >
         ::EllipticPreconditioner ParallelPreconditioner;
 
-    //Dune::SeqILU0<M,X,X>* ilu=new Dune::SeqILU0<M,X,X>(Ae, relax);
+    typedef typename CPRSelector<M,X,X,P>
-    typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
         ::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<M,X>( A_, param_.cpr_relax_, comm_ );
+                pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, comm, param_.cpr_relax_ );
             }
             else {
-                pre_ = ISTLUtility::createILUnPtr<M,X>( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm_ );
+                pre_ = ISTLUtility::createILUnPtr<M,X>( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_);
             }
         }
 

opm/autodiff/ISTLSolver.hpp

@@ -20,7 +20,6 @@
 #ifndef OPM_ISTLSOLVER_HEADER_INCLUDED
 #define OPM_ISTLSOLVER_HEADER_INCLUDED
 
-#include <opm/autodiff/AdditionalObjectDeleter.hpp>
 #include <opm/autodiff/BlackoilAmg.hpp>
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
@@ -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;

opm/autodiff/NewtonIterationBlackoilInterleaved.cpp

@@ -27,7 +27,6 @@
 #define FLOW_SUPPORT_AMG !defined(HAVE_UMFPACK)
 
 #include <opm/autodiff/DuneMatrix.hpp>
-#include <opm/autodiff/AdditionalObjectDeleter.hpp>
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
 #include <opm/autodiff/NewtonIterationUtilities.hpp>

opm/autodiff/ParallelOverlappingILU0.hpp

@@ -265,6 +265,26 @@ public:
         init( reinterpret_cast<const Matrix&>(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<class BlockType, class Alloc>
+    ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& 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<const Matrix&>(A), n );
+    }
+
     /*! \brief Constructor.
 
       Constructor gets all parameters to operate the prec.