From 221565f038705b4c4474df3a6ae66237d7d5bb60 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Wed, 28 Jan 2015 15:37:07 +0100
Subject: [PATCH 01/15] Enable the use of parallel dune-istl solvers.

As with opm-core we use boost::any to provide additional
information about a parallel run. It is used to set a
ParallelISTLInformation object and and fill it with the
information obtained from a parallel Cpgrid.

Note that the simulator currently compiles sucessfully. Still,
we have to test the runs and do debugging.
---
 CMakeLists_files.cmake                        |   2 +
 examples/sim_fibo_ad_cp.cpp                   |   8 +-
 opm/autodiff/CPRPreconditioner.hpp            | 243 ++++++++++++++++--
 .../ExtractParallelGridInformationToISTL.cpp  |  56 ++++
 .../ExtractParallelGridInformationToISTL.hpp  |  44 ++++
 opm/autodiff/NewtonIterationBlackoilCPR.cpp   |  61 +++--
 opm/autodiff/NewtonIterationBlackoilCPR.hpp   |  46 +++-
 .../NewtonIterationBlackoilSimple.cpp         |   9 +-
 .../NewtonIterationBlackoilSimple.hpp         |   6 +-
 9 files changed, 412 insertions(+), 63 deletions(-)
 create mode 100644 opm/autodiff/ExtractParallelGridInformationToISTL.cpp
 create mode 100644 opm/autodiff/ExtractParallelGridInformationToISTL.hpp

diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 6ee83b55e..9ebb3a6c4 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -28,6 +28,7 @@
 list (APPEND MAIN_SOURCE_FILES
 	opm/autodiff/BlackoilPropsAd.cpp
 	opm/autodiff/BlackoilPropsAdInterface.cpp
+	opm/autodiff/ExtractParallelGridInformationToISTL.cpp
 	opm/autodiff/NewtonIterationBlackoilCPR.cpp
 	opm/autodiff/NewtonIterationBlackoilSimple.cpp
 	opm/autodiff/GridHelpers.cpp
@@ -101,6 +102,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/CPRPreconditioner.hpp
 	opm/autodiff/fastSparseProduct.hpp
 	opm/autodiff/DuneMatrix.hpp
+  opm/autodiff/ExtractParallelGridInformationToISTL.hpp
 	opm/autodiff/GeoProps.hpp
 	opm/autodiff/GridHelpers.hpp
 	opm/autodiff/ImpesTPFAAD.hpp
diff --git a/examples/sim_fibo_ad_cp.cpp b/examples/sim_fibo_ad_cp.cpp
index c16f6a684..b309727b6 100644
--- a/examples/sim_fibo_ad_cp.cpp
+++ b/examples/sim_fibo_ad_cp.cpp
@@ -61,6 +61,7 @@
 #include <opm/core/linalg/LinearSolverFactory.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
+#include <opm/autodiff/ExtractParallelGridInformationToISTL.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
@@ -210,10 +211,13 @@ try
 
     // Solver for Newton iterations.
     std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
+    
+    boost::any parallel_information;
+    Opm::extractParallelGridInformationToISTL(parallel_information, *grid);
     if (param.getDefault("use_cpr", true)) {
-        fis_solver.reset(new NewtonIterationBlackoilCPR(param));
+        fis_solver.reset(new NewtonIterationBlackoilCPR(param, parallel_information));
     } else {
-        fis_solver.reset(new NewtonIterationBlackoilSimple(param));
+        fis_solver.reset(new NewtonIterationBlackoilSimple(param, parallel_information));
     }
 
     // Write parameters used for later reference.
diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp
index f7207c8a6..8b943a203 100644
--- a/opm/autodiff/CPRPreconditioner.hpp
+++ b/opm/autodiff/CPRPreconditioner.hpp
@@ -22,6 +22,7 @@
 #define OPM_CPRPRECONDITIONER_HEADER_INCLUDED
 
 #include <memory>
+#include <type_traits>
 
 #include <opm/core/utility/platform_dependent/disable_warnings.h>
 
@@ -43,6 +44,174 @@
 
 namespace Opm
 {
+namespace
+{
+//! \brief A custom deleter for the parallel preconditioners.
+//!
+//! In dune-istl they hold a reference to the sequential preconditioner.
+//! In CPRPreconditioner we use unique_ptr for the memory management.
+//! Ergo we need to construct the sequential preconditioner with new and
+//! make sure that it gets delete together with the enclosing parallel
+//! preconditioner. Therefore this delete store a pointer to it and deletes
+//! it during destruction.
+template<class PREC>
+class ParallelPreconditionerDeleter
+{
+public:
+    ParallelPreconditionerDeleter()
+        : ilu_()
+    {}
+    ParallelPreconditionerDeleter(PREC& ilu)
+    : ilu_(&ilu){}
+    template<class T>
+    void operator()(T* pt)
+    {
+        if(ilu_)
+        {
+            delete ilu_;
+        }
+        delete pt;
+    }
+private:
+    PREC* ilu_;
+};
+
+template<class M, class X, class Y, class P>
+struct CPRSelector
+{
+    typedef Dune::Amg::SequentialInformation ParallelInformation;
+    typedef Dune::MatrixAdapter<M, X, Y> Operator;
+    typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
+    typedef std::unique_ptr<EllipticPreconditioner> EllipticPreconditionerPointer;
+    
+    static Operator* makeOperator(const M& m, const P&)
+    {
+        return new Operator(m);
+    }
+};
+
+template<class M, class X, class Y, class I1, class I2>
+struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
+{
+    
+    typedef Dune::OwnerOverlapCopyCommunication<I1,I2> ParallelInformation;
+    typedef Dune::OverlappingSchwarzOperator<M,X,X,ParallelInformation> Operator;
+    typedef Dune::BlockPreconditioner<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
+    EllipticPreconditioner;
+    typedef std::unique_ptr<EllipticPreconditioner,
+                            ParallelPreconditionerDeleter<Dune::SeqILU0<M,X,X> > >
+    EllipticPreconditionerPointer;
+    
+    static Operator* makeOperator(const M& m, const ParallelInformation& p)
+    {
+        return new Operator(m, p);
+    }    
+};
+
+//! \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<class ILU, class I1, class I2>
+ParallelPreconditionerDeleter<ILU>
+createParallelDeleter(ILU& ilu, const Dune::OwnerOverlapCopyCommunication<I1,I2>& p)
+    {
+        (void) p;
+        return ParallelPreconditionerDeleter<ILU>(ilu);
+    }
+
+//! \brief Creates and initializes a unique pointer to an ILU0 preconditioner.
+template<class M, class X>
+std::shared_ptr<Dune::SeqILU0<M,X,X> >
+createILU0Ptr(const M& A_, double relax_, const Dune::Amg::SequentialInformation&)
+{
+    return std::shared_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>( A_, relax_) );
+}
+//! \brief Creates and initializes a shared pointer to an ILUn preconditioner.
+template<class M, class X>
+std::shared_ptr<Dune::SeqILUn<M,X,X> >
+createILUnPtr(const M& A_, int ilu_n, double relax_, const Dune::Amg::SequentialInformation&)
+{
+    return std::shared_ptr<Dune::SeqILUn<M,X,X> >(new Dune::SeqILUn<M,X,X>( A_, ilu_n, relax_) );
+}
+
+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.
+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 = new Dune::SeqILU0<M,X,X>( A_, relax_);
+    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 ILU0 preconditioner.
+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));
+}
+
+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_));
+}
+
+template<class M, class X=typename M::range_type, class I1, class I2>
+typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
+::EllipticPreconditionerPointer
+createEllipticPreconditionerPointer(const M& Ae, double relax_,
+                                    const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
+{
+    typedef Dune::BlockPreconditioner<X, X,
+                                      Dune::OwnerOverlapCopyCommunication<I1,I2>,
+                                      Dune::SeqILU0<M,X,X> >
+    ParallelPreconditioner;
+
+    Dune::SeqILU0<M,X,X>* ilu=new Dune::SeqILU0<M,X,X>(Ae, relax_);
+    typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
+        ::EllipticPreconditionerPointer EllipticPreconditionerPointer;
+    return EllipticPreconditionerPointer(new ParallelPreconditioner(*ilu, comm),
+                                         createParallelDeleter(*ilu, comm));
+}
+} // end namespace
+
+
     /*!
       \brief Sequential CPR preconditioner.
 
@@ -53,13 +222,16 @@ namespace Opm
       \tparam X Type of the update
       \tparam Y Type of the defect
     */
-    template<class M, class X, class Y>
+    template<class M, class X, class Y,
+             class P=Dune::Amg::SequentialInformation>
     class CPRPreconditioner : public Dune::Preconditioner<X,Y>
     {
         // prohibit copying for now
         CPRPreconditioner( const CPRPreconditioner& );
 
     public:
+        //! \brief The type describing the parallel information
+        typedef P ParallelInformation;
         //! \brief The matrix type the preconditioner is for.
         typedef typename Dune::remove_const<M>::type matrix_type;
         //! \brief The domain type of the preconditioner.
@@ -72,21 +244,28 @@ namespace Opm
         // define the category
         enum {
             //! \brief The category the preconditioner is part of.
-            category = Dune::SolverCategory::sequential
+            category = std::is_same<P,Dune::Amg::SequentialInformation>::value?
+            Dune::SolverCategory::sequential:Dune::SolverCategory::overlapping
         };
 
         //! \brief Elliptic Operator
-        typedef Dune::MatrixAdapter<M,X,X> Operator;
+        typedef typename CPRSelector<M,X,X,P>::Operator Operator;
 
         //! \brief preconditioner for the whole system (here either ILU(0) or ILU(n)
         typedef Dune::Preconditioner<X,X> WholeSystemPreconditioner;
 
-        //! \brief ilu-0 preconditioner for the elliptic system
-        typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
+        //! \brief the ilu-0 preconditioner used the for the elliptic system
+        typedef typename CPRSelector<M,X,X,P>::EllipticPreconditioner
+        EllipticPreconditioner;
+
+        //! \brief type of the unique pointer to the ilu-0 preconditioner
+        //! used the for the elliptic system
+        typedef typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
+        EllipticPreconditionerPointer;
 
         //! \brief amg preconditioner for the elliptic system
         typedef EllipticPreconditioner Smoother;
-        typedef Dune::Amg::AMG<Operator, X, Smoother> AMG;
+        typedef Dune::Amg::AMG<Operator, X, Smoother, P> AMG;
 
         /*! \brief Constructor.
 
@@ -96,32 +275,36 @@ namespace Opm
           \param relax   The ILU0 relaxation factor.
           \param useAMG  if true, AMG is used as a preconditioner for the elliptic sub-system, otherwise ilu-0 (default)
           \param useBiCG if true, BiCG solver is used (default), otherwise CG solver
+          \param paralleInformation The information about the parallelization, if this is a 
+                                    paralle run
         */
         CPRPreconditioner (const M& A, const M& Ae, const field_type relax,
                            const unsigned int ilu_n,
                            const bool useAMG,
-                           const bool useBiCG )
+                           const bool useBiCG,
+                           const ParallelInformation& comm=ParallelInformation())
             : A_(A),
               Ae_(Ae),
               de_( Ae_.N() ),
               ve_( Ae_.M() ),
               dmodified_( A_.N() ),
-              opAe_( Ae_ ),
+              opAe_(CPRSelector<M,X,Y,P>::makeOperator(Ae_, comm)),
               precond_(), // ilu0 preconditioner for elliptic system
               amg_(),     // amg  preconditioner for elliptic system
               pre_(), // copy A will be made be the preconditioner
               vilu_( A_.N() ),
               relax_(relax),
-              use_bicg_solver_( useBiCG )
+              use_bicg_solver_( useBiCG ),
+              comm_(comm)
         {
             // create appropriate preconditioner for elliptic system
-            createPreconditioner( useAMG );
+            createPreconditioner( useAMG, comm );
 
             if( ilu_n == 0 ) {
-                pre_.reset( new Dune::SeqILU0<M,X,X>( A_, relax_) );
+                pre_ = createILU0Ptr<M,X>( A_, relax_, comm );
             }
             else {
-                pre_.reset( new Dune::SeqILUn<M,X,X>( A_, ilu_n, relax_) );
+                pre_ = createILUnPtr<M,X>( A_, ilu_n, relax_, comm );
             }
         }
 
@@ -193,17 +376,22 @@ namespace Opm
             // operator result containing iterations etc.
             Dune::InverseOperatorResult result;
 
-            // sequential scalar product
-            Dune::SeqScalarProduct<X> sp;
+            // the scalar product chooser
+            typedef Dune::ScalarProductChooser<X,ParallelInformation,category>
+                ScalarProductChooser;
+            // the scalar product.
+            std::unique_ptr<typename ScalarProductChooser::ScalarProduct> 
+                sp(ScalarProductChooser::construct(comm_));
+
             if( amg_ )
             {
                 // Solve system with AMG
                 if( use_bicg_solver_ ) {
-                    Dune::BiCGSTABSolver<X> linsolve(opAe_, sp, (*amg_), tolerance, maxit, verbosity);
+                    Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
                     linsolve.apply(x, de, result);
                 }
                 else {
-                    Dune::CGSolver<X> linsolve(opAe_, sp, (*amg_), tolerance, maxit, verbosity);
+                    Dune::CGSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
                     linsolve.apply(x, de, result);
                 }
             }
@@ -212,11 +400,11 @@ namespace Opm
                 assert( precond_ );
                 // Solve system with ILU-0
                 if( use_bicg_solver_ ) {
-                    Dune::BiCGSTABSolver<X> linsolve(opAe_, sp, (*precond_), tolerance, maxit, verbosity);
+                    Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
                     linsolve.apply(x, de, result);
                 }
                 else {
-                    Dune::CGSolver<X> linsolve(opAe_, sp, (*precond_), tolerance, maxit, verbosity);
+                    Dune::CGSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
                     linsolve.apply(x, de, result);
                 }
 
@@ -236,15 +424,20 @@ namespace Opm
         Y de_, ve_, dmodified_;
 
         //! \brief elliptic operator
-        Operator opAe_;
+        std::unique_ptr<Operator> opAe_;
 
         //! \brief ILU0 preconditioner for the elliptic system
-        std::unique_ptr< EllipticPreconditioner > precond_;
+        EllipticPreconditionerPointer precond_;
         //! \brief AMG preconditioner with ILU0 smoother
         std::unique_ptr< AMG > amg_;
 
         //! \brief The preconditioner for the whole system
-        std::unique_ptr< WholeSystemPreconditioner > pre_;
+        //!
+        //! We have to use a shared_ptr instead of a unique_ptr
+        //! as we need to use a custom allocator based on dynamic
+        //! information. But for unique_ptr the type of this deleter
+        //! has to be available at coompile time.
+        std::shared_ptr< WholeSystemPreconditioner > pre_;
 
         //! \brief temporary variables for ILU solve
         Y vilu_;
@@ -255,8 +448,10 @@ namespace Opm
         //! \brief true if ISTL BiCGSTABSolver is used, otherwise ISTL CGSolver is used
         const bool use_bicg_solver_;
 
+        //! \brief The information about the parallelization
+        const P& comm_;
      protected:
-        void createPreconditioner( const bool amg )
+        void createPreconditioner( const bool amg, const P& comm )
         {
             if( amg )
             {
@@ -275,10 +470,10 @@ namespace Opm
               criterion.setAlpha(.67);
               criterion.setBeta(1.0e-6);
               criterion.setMaxLevel(10);
-              amg_ = std::unique_ptr< AMG > (new AMG(opAe_, criterion, smootherArgs));
+              amg_ = std::unique_ptr< AMG > (new AMG(*opAe_, criterion, smootherArgs));
             }
             else
-              precond_ = std::unique_ptr< EllipticPreconditioner > (new EllipticPreconditioner( Ae_, relax_ ));
+                precond_ = createEllipticPreconditionerPointer<M,X>( Ae_, relax_, comm);
        }
     };
 
diff --git a/opm/autodiff/ExtractParallelGridInformationToISTL.cpp b/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
new file mode 100644
index 000000000..42f30a220
--- /dev/null
+++ b/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
@@ -0,0 +1,56 @@
+/*
+  Copyright 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
+
+  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/>.
+*/
+#include <config.h>
+#include "ExtractParallelGridInformationToISTL.hpp"
+#include <opm/core/linalg/ParallelIstlInformation.hpp>
+#include <dune/common/version.hh>
+
+ namespace Opm
+{
+#if defined(HAVE_DUNE_CORNERPOINT) && DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
+#if defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL)
+// Extracts the information about the data decomposition from the grid for dune-istl
+struct NullDeleter
+{
+    void operator()(void*)
+    {}
+};
+
+void extractParallelGridInformationToISTL(boost::any& anyComm, const Dune::CpGrid& grid)
+{
+    if(grid.comm().size()>1)
+    {
+        // this is a parallel run with distributed data.
+        NullDeleter null_delete;
+        Dune::CpGrid& mgrid=const_cast<Dune::CpGrid&>(grid);
+        Dune::CpGrid::ParallelIndexSet* idx=&(mgrid.getCellIndexSet());
+        Dune::CpGrid::RemoteIndices* ridx=&(mgrid.getCellRemoteIndices());
+        anyComm=boost::any(Opm::ParallelISTLInformation(std::shared_ptr<Dune::CpGrid::ParallelIndexSet>(idx, null_delete),
+                                                        std::shared_ptr<Dune::CpGrid::RemoteIndices>(ridx, null_delete),
+                                                        grid.comm()));
+    }
+}
+#else
+// Missing support for MPI or dune-istl -> do nothing.
+void extractParallelGridInformationToISTL(boost::any&, const Dune::CpGrid&)
+{}
+#endif //defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL)
+#endif //defined(HAVE_DUNE_CORNERPOINT)
+} // end namespace Opm
diff --git a/opm/autodiff/ExtractParallelGridInformationToISTL.hpp b/opm/autodiff/ExtractParallelGridInformationToISTL.hpp
new file mode 100644
index 000000000..4f4b08e40
--- /dev/null
+++ b/opm/autodiff/ExtractParallelGridInformationToISTL.hpp
@@ -0,0 +1,44 @@
+/*
+  Copyright 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
+
+  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_EXTRACTPARALLELGRIDINFORMATIONTOISTL_HEADER_INCLUDED
+#define OPM_EXTRACTPARALLELGRIDINFORMATIONTOISTL_HEADER_INCLUDED
+#ifdef HAVE_DUNE_CORNERPOINT
+
+#include<dune/grid/CpGrid.hpp>
+#include<boost/any.hpp>
+
+namespace Opm
+{
+
+/// \brief Extracts the information about the data decomposition from the grid for dune-istl
+/// 
+/// In the case that grid is a parallel grid this method will query it to get the information
+/// about the data decompoisition and convert it to the format expected by the linear algebra
+/// of dune-istl.
+/// \warn if there is no support for dune-istl and MPI then this functio does not do anything.
+/// \param anyComm The handle to to store the information in. If grid is a parallel grid
+/// then this will ecapsulate an instance of ParallelISTLInformation.
+/// \param grid The grid to inspect.
+
+void extractParallelGridInformationToISTL(boost::any& anyComm, const Dune::CpGrid& grid);
+
+} // end namespace Opm
+#endif //defined(HAVE_DUNE_CORNERPOINT)
+#endif // OPM_EXTRACTPARALLELGRIDINFORMATIONTOISTL_HEADER_INCLUDED
diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
index 0eb295e8b..80933b2b7 100644
--- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp
+++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
@@ -22,17 +22,14 @@
 #include <opm/autodiff/DuneMatrix.hpp>
 
 #include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
-#include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/Units.hpp>
 #include <opm/core/linalg/LinearSolverFactory.hpp>
-
+#include <opm/core/linalg/ParallelIstlInformation.hpp>
 #include <opm/core/utility/platform_dependent/disable_warnings.h>
 
-#include <dune/istl/bvector.hh>
 // #include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/operators.hh>
 #include <dune/istl/io.hh>
 #include <dune/istl/owneroverlapcopy.hh>
 #include <dune/istl/preconditioners.hh>
@@ -57,11 +54,8 @@ namespace Opm
     typedef AutoDiffBlock<double> ADB;
     typedef ADB::V V;
     typedef ADB::M M;
-
-    typedef Dune::FieldVector<double, 1   > VectorBlockType;
     typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
     typedef Dune::BCRSMatrix <MatrixBlockType>        Mat;
-    typedef Dune::BlockVector<VectorBlockType>        Vector;
 
 
     namespace {
@@ -114,8 +108,9 @@ namespace Opm
 
 
     /// Construct a system solver.
-    NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param)
-        : iterations_( 0 )
+    NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param,
+                                                           const boost::any& parallelInformation)
+        : iterations_( 0 ), parallelInformation_(parallelInformation)
     {
         cpr_relax_        = param.getDefault("cpr_relax", 1.0);
         cpr_ilu_n_        = param.getDefault("cpr_ilu_n", 0);
@@ -186,34 +181,38 @@ namespace Opm
 
         // Create ISTL matrix for elliptic part.
         DuneMatrix istlAe( A.topLeftCorner(nc, nc) );
-
-        // Construct operator, scalar product and vectors needed.
-        typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
-        Operator opA(istlA);
-        Dune::SeqScalarProduct<Vector> sp;
+        
         // Right hand side.
-        Vector istlb(opA.getmat().N());
+        Vector istlb(istlA.N());
         std::copy_n(b.data(), istlb.size(), istlb.begin());
         // System solution
-        Vector x(opA.getmat().M());
+        Vector x(istlA.M());
         x = 0.0;
 
-        // Construct preconditioner.
-        // typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
-        typedef Opm::CPRPreconditioner<Mat,Vector,Vector> Preconditioner;
-        Preconditioner precond(istlA, istlAe, cpr_relax_, cpr_ilu_n_, cpr_use_amg_, cpr_use_bicgstab_);
-
-        // Construct linear solver.
-        const double tolerance = 1e-3;
-        const int maxit = 150;
-        const int verbosity = 0;
-        const int restart = 40;
-        Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond, tolerance, restart, maxit, verbosity);
-
-        // Solve system.
         Dune::InverseOperatorResult result;
-        linsolve.apply(x, istlb, result);
-
+#if HAVE_MPI
+        if(parallelInformation_.type()==typeid(ParallelISTLInformation))
+        {
+            typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
+            const ParallelISTLInformation& info =
+                boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
+            Comm istlComm(info.communicator());
+            info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices());
+            // Construct operator, scalar product and vectors needed.
+            typedef Dune::OverlappingSchwarzOperator<Mat,Vector,Vector,Comm> Operator;
+            Operator opA(istlA, istlComm);
+            constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, istlAe, x, istlb, istlComm, result);
+        }
+        else
+#endif
+        {
+            // Construct operator, scalar product and vectors needed.
+            typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
+            Operator opA(istlA);
+            Dune::Amg::SequentialInformation info;
+            constructPreconditionerAndSolve(opA, istlAe, x, istlb, info, result);
+        }
+            
         // store number of iterations
         iterations_ = result.iterations;
 
diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.hpp b/opm/autodiff/NewtonIterationBlackoilCPR.hpp
index ae73bed46..58a19989d 100644
--- a/opm/autodiff/NewtonIterationBlackoilCPR.hpp
+++ b/opm/autodiff/NewtonIterationBlackoilCPR.hpp
@@ -20,10 +20,14 @@
 #ifndef OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED
 #define OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED
 
-
+#include <opm/autodiff/DuneMatrix.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
+#include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <dune/istl/scalarproducts.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/bvector.hh>
 #include <memory>
 
 namespace Opm
@@ -37,7 +41,12 @@ namespace Opm
     /// in Fully Implicit Reservoir Simulations" by Gries et al (SPE 163608).
     class NewtonIterationBlackoilCPR : public NewtonIterationBlackoilInterface
     {
+        typedef Dune::FieldVector<double, 1   > VectorBlockType;
+        typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
+        typedef Dune::BCRSMatrix <MatrixBlockType>        Mat;
+        typedef Dune::BlockVector<VectorBlockType>        Vector;
     public:
+        
         /// Construct a system solver.
         /// \param[in] param   parameters controlling the behaviour of
         ///                    the preconditioning and choice of
@@ -47,7 +56,10 @@ namespace Opm
         ///                        cpr_ilu_n        (default 0) use ILU(n) for preconditioning of the linear system
         ///                        cpr_use_amg      (default false) if true, use AMG preconditioner for elliptic part
         ///                        cpr_use_bicgstab (default true)  if true, use BiCGStab (else use CG) for elliptic part
-        NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param);
+        /// \param[in] parallelInformation In the case of a parallel run
+        ///                               with dune-istl the information about the parallelization.
+        NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param,
+                                   const boost::any& parallelInformation=boost::any());
 
         /// Solve the system of linear equations Ax = b, with A being the
         /// combined derivative matrix of the residual and b
@@ -58,12 +70,42 @@ namespace Opm
 
         /// \copydoc NewtonIterationBlackoilInterface::iterations
         virtual int iterations () const { return iterations_; }
+
+        /// \brief construct the CPR preconditioner and the solver.
+        /// \tparam P The type of the parallel information.
+        /// \param parallelInformation the information about the parallelization.
+        template<int category=Dune::SolverCategory::sequential, class O, class P>   
+        void constructPreconditionerAndSolve(O& opA, DuneMatrix istlAe, 
+                                             Vector& x, Vector& istlb,
+                                             const P& parallelInformation,
+                                             Dune::InverseOperatorResult& result) const
+        {
+            typedef Dune::ScalarProductChooser<Vector,P,category> ScalarProductChooser;
+            std::unique_ptr<typename ScalarProductChooser::ScalarProduct> 
+                sp(ScalarProductChooser::construct(parallelInformation));
+            // Construct preconditioner.
+            // typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
+            typedef Opm::CPRPreconditioner<Mat,Vector,Vector,P> Preconditioner;
+            Preconditioner precond(opA.getmat(), istlAe, cpr_relax_, cpr_ilu_n_, cpr_use_amg_, cpr_use_bicgstab_, parallelInformation);
+
+            // Construct linear solver.
+            const double tolerance = 1e-3;
+            const int maxit = 150;
+            const int verbosity = 0;
+            const int restart = 40;
+            Dune::RestartedGMResSolver<Vector> linsolve(opA, *sp, precond, tolerance, restart, maxit, verbosity);
+            
+            // Solve system.
+            linsolve.apply(x, istlb, result);
+        }
+
     private:
         mutable int iterations_;
         double cpr_relax_;
         unsigned int cpr_ilu_n_;
         bool cpr_use_amg_;
         bool cpr_use_bicgstab_;
+        const boost::any& parallelInformation_;
     };
 
 } // namespace Opm
diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.cpp b/opm/autodiff/NewtonIterationBlackoilSimple.cpp
index 228225b4c..b976680b4 100644
--- a/opm/autodiff/NewtonIterationBlackoilSimple.cpp
+++ b/opm/autodiff/NewtonIterationBlackoilSimple.cpp
@@ -29,8 +29,11 @@ namespace Opm
 
     /// Construct a system solver.
     /// \param[in] linsolver   linear solver to use
-    NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param)
-        : iterations_( 0 )
+    /// \param[in] parallelInformation In the case of a parallel run
+    ///  with dune-istl the information about the parallelization.
+    NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param,
+                                                                 const boost::any& parallelInformation)
+        : iterations_( 0 ), parallelInformation_(parallelInformation)
     {
         linsolver_.reset(new LinearSolverFactory(param));
     }
@@ -58,7 +61,7 @@ namespace Opm
         Opm::LinearSolverInterface::LinearSolverReport rep
             = linsolver_->solve(matr.rows(), matr.nonZeros(),
                                 matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
-                                total_residual.value().data(), dx.data());
+                                total_residual.value().data(), dx.data(), parallelInformation_);
 
         // store iterations
         iterations_ = rep.iterations;
diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.hpp b/opm/autodiff/NewtonIterationBlackoilSimple.hpp
index 6bdaf7b07..20cd484c4 100644
--- a/opm/autodiff/NewtonIterationBlackoilSimple.hpp
+++ b/opm/autodiff/NewtonIterationBlackoilSimple.hpp
@@ -40,7 +40,10 @@ namespace Opm
         /// Construct a system solver.
         /// \param[in] param   parameters controlling the behaviour and
         ///                    choice of linear solver.
-        NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param);
+        /// \param[in] parallelInformation In the case of a parallel run
+        ///                    with dune-istl the information about the parallelization.
+        NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param,
+                                      const boost::any& parallelInformation=boost::any());
 
         /// Solve the system of linear equations Ax = b, with A being the
         /// combined derivative matrix of the residual and b
@@ -55,6 +58,7 @@ namespace Opm
     private:
         std::unique_ptr<LinearSolverInterface> linsolver_;
         mutable int iterations_;
+        const boost::any& parallelInformation_;
     };
 
 } // namespace Opm

From 4527ce8ffd61d79167c4498f184cdb76a92743ac Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Thu, 29 Jan 2015 11:30:03 +0100
Subject: [PATCH 02/15] Add access to the underlying information about the
 parallelization.

We need it serveral places and all of them seem to have access to
NewtonIterationBlackoilInterface. This makes it natural to give access
to it and prevent users from having to forward it manually at several
places in the simulator driver.
---
 opm/autodiff/NewtonIterationBlackoilCPR.cpp       | 4 ++++
 opm/autodiff/NewtonIterationBlackoilCPR.hpp       | 6 +++++-
 opm/autodiff/NewtonIterationBlackoilInterface.hpp | 5 +++++
 opm/autodiff/NewtonIterationBlackoilSimple.cpp    | 4 ++++
 opm/autodiff/NewtonIterationBlackoilSimple.hpp    | 3 +++
 5 files changed, 21 insertions(+), 1 deletion(-)

diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
index 80933b2b7..de76309a0 100644
--- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp
+++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
@@ -234,6 +234,10 @@ namespace Opm
         return dx;
     }
 
+    const boost::any& NewtonIterationBlackoilCPR::parallelInformation() const
+    {
+        return parallelInformation_;
+    }
 
 
 
diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.hpp b/opm/autodiff/NewtonIterationBlackoilCPR.hpp
index 58a19989d..e97a56408 100644
--- a/opm/autodiff/NewtonIterationBlackoilCPR.hpp
+++ b/opm/autodiff/NewtonIterationBlackoilCPR.hpp
@@ -71,6 +71,11 @@ namespace Opm
         /// \copydoc NewtonIterationBlackoilInterface::iterations
         virtual int iterations () const { return iterations_; }
 
+        /// \copydoc NewtonIterationBlackoilInterface::parallelInformation
+        virtual const boost::any& parallelInformation() const;
+
+    private:
+
         /// \brief construct the CPR preconditioner and the solver.
         /// \tparam P The type of the parallel information.
         /// \param parallelInformation the information about the parallelization.
@@ -99,7 +104,6 @@ namespace Opm
             linsolve.apply(x, istlb, result);
         }
 
-    private:
         mutable int iterations_;
         double cpr_relax_;
         unsigned int cpr_ilu_n_;
diff --git a/opm/autodiff/NewtonIterationBlackoilInterface.hpp b/opm/autodiff/NewtonIterationBlackoilInterface.hpp
index 2e26db07e..ffe886153 100644
--- a/opm/autodiff/NewtonIterationBlackoilInterface.hpp
+++ b/opm/autodiff/NewtonIterationBlackoilInterface.hpp
@@ -23,6 +23,7 @@
 
 #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
 
+#include <boost/any.hpp>
 namespace Opm
 {
 
@@ -42,6 +43,10 @@ namespace Opm
 
         /// \return number of linear iterations used during last call of computeNewtonIncrement
         virtual int iterations () const = 0;
+
+
+        /// \brief Get the information about the parallelization of the grid.
+        virtual const boost::any& parallelInformation() const = 0;
     };
 
 } // namespace Opm
diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.cpp b/opm/autodiff/NewtonIterationBlackoilSimple.cpp
index b976680b4..3817d94b9 100644
--- a/opm/autodiff/NewtonIterationBlackoilSimple.cpp
+++ b/opm/autodiff/NewtonIterationBlackoilSimple.cpp
@@ -74,5 +74,9 @@ namespace Opm
         return dx;
     }
 
+    const boost::any& NewtonIterationBlackoilSimple::parallelInformation() const
+    {
+        return parallelInformation_;
+    }
 } // namespace Opm
 
diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.hpp b/opm/autodiff/NewtonIterationBlackoilSimple.hpp
index 20cd484c4..5f530f10c 100644
--- a/opm/autodiff/NewtonIterationBlackoilSimple.hpp
+++ b/opm/autodiff/NewtonIterationBlackoilSimple.hpp
@@ -55,6 +55,9 @@ namespace Opm
         /// \copydoc NewtonIterationBlackoilInterface::iterations
         virtual int iterations () const { return iterations_; }
 
+        /// \copydoc NewtonIterationBlackoilInterface::parallelInformation
+        virtual const boost::any& parallelInformation() const;
+
     private:
         std::unique_ptr<LinearSolverInterface> linsolver_;
         mutable int iterations_;

From ff9b8d790d6415827e4d880c53f2a14df8b72b45 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Thu, 29 Jan 2015 11:32:39 +0100
Subject: [PATCH 03/15] Added a parallel version for computing the global
 reductions.

---
 .../FullyImplicitBlackoilSolver_impl.hpp      | 48 +++++++++++++++----
 1 file changed, 40 insertions(+), 8 deletions(-)

diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index dc23a1fb3..a04ad458a 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -29,6 +29,7 @@
 
 #include <opm/core/grid.h>
 #include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <opm/core/linalg/ParallelIstlInformation.hpp>
 #include <opm/core/props/rock/RockCompressibility.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
@@ -1898,16 +1899,47 @@ namespace {
                                                          int nc) const
     {
         // Do the global reductions
-        for ( int idx=0; idx<MaxNumPhases; ++idx )
+        #if HAVE_MPI
+        if(linsolver_.parallelInformation().type()==typeid(ParallelISTLInformation))
         {
-            if (active_[idx]) {
-                B_avg[idx] = B.col(idx).sum()/nc;
-                maxCoeff[idx]=tempV.col(idx).maxCoeff();
-                R_sum[idx] = R.col(idx).sum();
-            }
-            else
+            const ParallelISTLInformation& info =
+                boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
+            for ( int idx=0; idx<MaxNumPhases; ++idx )
             {
-                R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.;
+                if (active_[idx]) {
+                    auto values     = std::tuple<double,double,double>(0.,0.,0.);
+                    auto containers = std::make_tuple(B.col(idx),
+                                                      tempV.col(idx),
+                                                      R.col(idx));
+                    auto operators  = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
+                                                      Opm::Reduction::makeGlobalMaxFunctor<double>(),
+                                                      Opm::Reduction::makeGlobalSumFunctor<double>());
+                    nc=info.communicator().sum(nc);
+                    info.computeReduction(containers,operators,values);
+                    B_avg[idx]    = std::get<0>(values)/nc;
+                    maxCoeff[idx] = std::get<1>(values);
+                    R_sum[idx]    = std::get<2>(values);
+                }
+                else
+                {
+                    R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.;
+                }
+            }
+        }
+        else
+#endif
+        {
+            for ( int idx=0; idx<MaxNumPhases; ++idx )
+            {
+                if (active_[idx]) {
+                    B_avg[idx] = B.col(idx).sum()/nc;
+                    maxCoeff[idx]=tempV.col(idx).maxCoeff();
+                    R_sum[idx] = R.col(idx).sum();
+                }
+                else
+                {
+                    R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.;
+                }
             }
         }
         // Compute total pore volume

From a7f1614f63bdcbf738c2267d4f90f5040ebfb5d4 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Tue, 3 Feb 2015 11:31:57 +0100
Subject: [PATCH 04/15] Adds copyright declarations and more documentation.

---
 opm/autodiff/CPRPreconditioner.hpp            | 96 ++++++++++++++-----
 .../ExtractParallelGridInformationToISTL.hpp  |  2 +-
 .../FullyImplicitBlackoilSolver_impl.hpp      |  2 +
 opm/autodiff/NewtonIterationBlackoilCPR.cpp   |  2 +
 4 files changed, 76 insertions(+), 26 deletions(-)

diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp
index 8b943a203..af35d5e82 100644
--- a/opm/autodiff/CPRPreconditioner.hpp
+++ b/opm/autodiff/CPRPreconditioner.hpp
@@ -1,6 +1,8 @@
 /*
   Copyright 2014 SINTEF ICT, Applied Mathematics.
   Copyright 2014 IRIS AS.
+  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -75,37 +77,57 @@ public:
 private:
     PREC* ilu_;
 };
-
+///
+/// \brief A traits class for selecting the types of the preconditioner.
+///
+/// \tparam M The type of the matrix.
+/// \tparam X The type of the domain of the linear problem.
+/// \tparam Y The type of the range of the linear problem.
+/// \tparam P The type of the parallel information.
+////
 template<class M, class X, class Y, class P>
 struct CPRSelector
 {
+    /// \brief The information about the parallelization and communication
     typedef Dune::Amg::SequentialInformation ParallelInformation;
+    /// \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 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);
     }
 };
 
+/// \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;
+    /// \brief The type of the unique pointer to the preconditioner of the elliptic part.
     typedef std::unique_ptr<EllipticPreconditioner,
                             ParallelPreconditionerDeleter<Dune::SeqILU0<M,X,X> > >
     EllipticPreconditionerPointer;
-    
+
+    /// \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& p)
     {
         return new Operator(m, p);
-    }    
+    }
 };
 
 //! \brief Creates the deleter needed for the parallel ILU preconditioners.
@@ -122,19 +144,24 @@ createParallelDeleter(ILU& ilu, const Dune::OwnerOverlapCopyCommunication<I1,I2>
         return ParallelPreconditionerDeleter<ILU>(ilu);
     }
 
-//! \brief Creates and initializes a unique pointer to an ILU0 preconditioner.
+//! \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>
 std::shared_ptr<Dune::SeqILU0<M,X,X> >
-createILU0Ptr(const M& A_, double relax_, const Dune::Amg::SequentialInformation&)
+createILU0Ptr(const M& A, double relax, const Dune::Amg::SequentialInformation&)
 {
-    return std::shared_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>( A_, relax_) );
+    return std::shared_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>( A, 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>
 std::shared_ptr<Dune::SeqILUn<M,X,X> >
-createILUnPtr(const M& A_, int ilu_n, double relax_, const Dune::Amg::SequentialInformation&)
+createILUnPtr(const M& A, int ilu_n, double relax, const Dune::Amg::SequentialInformation&)
 {
-    return std::shared_ptr<Dune::SeqILUn<M,X,X> >(new Dune::SeqILUn<M,X,X>( A_, ilu_n, relax_) );
+    return std::shared_ptr<Dune::SeqILUn<M,X,X> >(new Dune::SeqILUn<M,X,X>( A, ilu_n, relax) );
 }
 
 template<class ILU, class I1, class I2>
@@ -151,9 +178,12 @@ struct SelectParallelILUSharedPtr
 };
 
 //! \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_, 
+createILU0Ptr(const M& A, double relax,
               const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
 {
     typedef Dune::BlockPreconditioner<
@@ -162,40 +192,52 @@ createILU0Ptr(const M& A_, double relax_,
         Dune::OwnerOverlapCopyCommunication<I1,I2>,
         Dune::SeqILU0<M,X,X>
         > PointerType;
-    Dune::SeqILU0<M,X,X>* ilu = new Dune::SeqILU0<M,X,X>( A_, relax_);
+    Dune::SeqILU0<M,X,X>* ilu = new Dune::SeqILU0<M,X,X>(A, relax);
     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 ILU0 preconditioner.
+
+//! \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_, 
+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>    
+        Dune::SeqILUn<M,X,X>
         > PointerType;
-    Dune::SeqILUn<M,X,X>* ilu = new Dune::SeqILUn<M,X,X>( A_, ilu_n, relax_);
-        
+    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));
 }
 
+/// \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_,
+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_));
+    return std::unique_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>(Ae, relax));
 }
 
+/// \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>
 typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
 ::EllipticPreconditionerPointer
-createEllipticPreconditionerPointer(const M& Ae, double relax_,
+createEllipticPreconditionerPointer(const M& Ae, double relax,
                                     const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
 {
     typedef Dune::BlockPreconditioner<X, X,
@@ -203,7 +245,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax_,
                                       Dune::SeqILU0<M,X,X> >
     ParallelPreconditioner;
 
-    Dune::SeqILU0<M,X,X>* ilu=new Dune::SeqILU0<M,X,X>(Ae, relax_);
+    Dune::SeqILU0<M,X,X>* ilu=new Dune::SeqILU0<M,X,X>(Ae, relax);
     typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
         ::EllipticPreconditionerPointer EllipticPreconditionerPointer;
     return EllipticPreconditionerPointer(new ParallelPreconditioner(*ilu, comm),
@@ -213,7 +255,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax_,
 
 
     /*!
-      \brief Sequential CPR preconditioner.
+      \brief CPR preconditioner.
 
       This is a two-stage preconditioner, combining an elliptic-type
       partial solution with ILU0 for the whole system.
@@ -221,6 +263,10 @@ createEllipticPreconditionerPointer(const M& Ae, double relax_,
       \tparam M The matrix type to operate on
       \tparam X Type of the update
       \tparam Y Type of the defect
+      \tparam P Type of the parallel information. If not provided
+                this will be Dune::Amg::SequentialInformation.
+                The preconditioner is parallel if this is
+                Dune::OwnerOverlapCopyCommunication<int,int>
     */
     template<class M, class X, class Y,
              class P=Dune::Amg::SequentialInformation>
@@ -275,8 +321,8 @@ createEllipticPreconditionerPointer(const M& Ae, double relax_,
           \param relax   The ILU0 relaxation factor.
           \param useAMG  if true, AMG is used as a preconditioner for the elliptic sub-system, otherwise ilu-0 (default)
           \param useBiCG if true, BiCG solver is used (default), otherwise CG solver
-          \param paralleInformation The information about the parallelization, if this is a 
-                                    paralle run
+          \param paralleInformation The information about the parallelization, if this is a
+                                    parallel run
         */
         CPRPreconditioner (const M& A, const M& Ae, const field_type relax,
                            const unsigned int ilu_n,
@@ -380,7 +426,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax_,
             typedef Dune::ScalarProductChooser<X,ParallelInformation,category>
                 ScalarProductChooser;
             // the scalar product.
-            std::unique_ptr<typename ScalarProductChooser::ScalarProduct> 
+            std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
                 sp(ScalarProductChooser::construct(comm_));
 
             if( amg_ )
diff --git a/opm/autodiff/ExtractParallelGridInformationToISTL.hpp b/opm/autodiff/ExtractParallelGridInformationToISTL.hpp
index 4f4b08e40..ab57d456d 100644
--- a/opm/autodiff/ExtractParallelGridInformationToISTL.hpp
+++ b/opm/autodiff/ExtractParallelGridInformationToISTL.hpp
@@ -28,7 +28,7 @@ namespace Opm
 {
 
 /// \brief Extracts the information about the data decomposition from the grid for dune-istl
-/// 
+///
 /// In the case that grid is a parallel grid this method will query it to get the information
 /// about the data decompoisition and convert it to the format expected by the linear algebra
 /// of dune-istl.
diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index a04ad458a..7c262ea85 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -1,5 +1,7 @@
 /*
   Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
 
   This file is part of the Open Porous Media project (OPM).
 
diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
index de76309a0..11107f5e1 100644
--- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp
+++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
@@ -1,5 +1,7 @@
 /*
   Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
 
   This file is part of the Open Porous Media project (OPM).
 

From b73f1c86fc36122fae22589271b44ec75620a8cc Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Tue, 3 Feb 2015 11:42:05 +0100
Subject: [PATCH 05/15] Fixed typos in documentatio of deleter.

---
 opm/autodiff/CPRPreconditioner.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp
index af35d5e82..00f63bdc4 100644
--- a/opm/autodiff/CPRPreconditioner.hpp
+++ b/opm/autodiff/CPRPreconditioner.hpp
@@ -53,8 +53,8 @@ namespace
 //! In dune-istl they hold a reference to the sequential preconditioner.
 //! In CPRPreconditioner we use unique_ptr for the memory management.
 //! Ergo we need to construct the sequential preconditioner with new and
-//! make sure that it gets delete together with the enclosing parallel
-//! preconditioner. Therefore this delete store a pointer to it and deletes
+//! 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.
 template<class PREC>
 class ParallelPreconditionerDeleter

From 70d5d1856087251321f52377ad82d75d647d824a Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Wed, 4 Feb 2015 22:38:10 +0100
Subject: [PATCH 06/15] Prevent copying of a matrix.

---
 opm/autodiff/NewtonIterationBlackoilCPR.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.hpp b/opm/autodiff/NewtonIterationBlackoilCPR.hpp
index e97a56408..ccd2caa43 100644
--- a/opm/autodiff/NewtonIterationBlackoilCPR.hpp
+++ b/opm/autodiff/NewtonIterationBlackoilCPR.hpp
@@ -80,7 +80,7 @@ namespace Opm
         /// \tparam P The type of the parallel information.
         /// \param parallelInformation the information about the parallelization.
         template<int category=Dune::SolverCategory::sequential, class O, class P>   
-        void constructPreconditionerAndSolve(O& opA, DuneMatrix istlAe, 
+        void constructPreconditionerAndSolve(O& opA, DuneMatrix& istlAe,
                                              Vector& x, Vector& istlb,
                                              const P& parallelInformation,
                                              Dune::InverseOperatorResult& result) const

From 11cb82e815233a94559cd0abe6747026e6193bce Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Thu, 5 Feb 2015 09:47:44 +0100
Subject: [PATCH 07/15] Correctly compute the global number of cells and pore
 volume.

---
 .../FullyImplicitBlackoilSolver_impl.hpp      | 19 ++++++++++++++-----
 1 file changed, 14 insertions(+), 5 deletions(-)

diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index 7c262ea85..8a2883699 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -1901,11 +1901,19 @@ namespace {
                                                          int nc) const
     {
         // Do the global reductions
-        #if HAVE_MPI
+#if HAVE_MPI
         if(linsolver_.parallelInformation().type()==typeid(ParallelISTLInformation))
         {
             const ParallelISTLInformation& info =
                 boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
+            // Compute the global number of cells and porevolume
+            std::vector<int> v(nc, 1);
+            auto nc_and_pv = std::tuple<int, double>(0,.0);
+            auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
+                                                        Opm::Reduction::makeGlobalSumFunctor<double>());
+            auto nc_and_pv_containers  = std::make_tuple(v, geo_.poreVolume());
+            info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
+            
             for ( int idx=0; idx<MaxNumPhases; ++idx )
             {
                 if (active_[idx]) {
@@ -1916,9 +1924,8 @@ namespace {
                     auto operators  = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
                                                       Opm::Reduction::makeGlobalMaxFunctor<double>(),
                                                       Opm::Reduction::makeGlobalSumFunctor<double>());
-                    nc=info.communicator().sum(nc);
                     info.computeReduction(containers,operators,values);
-                    B_avg[idx]    = std::get<0>(values)/nc;
+                    B_avg[idx]    = std::get<0>(values)/std::get<0>(nc_and_pv);
                     maxCoeff[idx] = std::get<1>(values);
                     R_sum[idx]    = std::get<2>(values);
                 }
@@ -1927,6 +1934,8 @@ namespace {
                     R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.;
                 }
             }
+            // Compute pore volume
+            return std::get<1>(nc_and_pv);
         }
         else
 #endif
@@ -1943,9 +1952,9 @@ namespace {
                     R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.;
                 }
             }
+            // Compute total pore volume
+            return geo_.poreVolume().sum();
         }
-        // Compute total pore volume
-        return geo_.poreVolume().sum();
     }
 
     template<class T>

From 39a6e1909977abb5316a55d8c6f14df864e26996 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Thu, 5 Feb 2015 10:08:47 +0100
Subject: [PATCH 08/15] Fixes compilation issues when no MPI is available.

---
 opm/autodiff/CPRPreconditioner.hpp | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp
index 00f63bdc4..1d078345a 100644
--- a/opm/autodiff/CPRPreconditioner.hpp
+++ b/opm/autodiff/CPRPreconditioner.hpp
@@ -105,6 +105,7 @@ struct CPRSelector
     }
 };
 
+#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> >
@@ -143,6 +144,7 @@ createParallelDeleter(ILU& ilu, const Dune::OwnerOverlapCopyCommunication<I1,I2>
         (void) p;
         return ParallelPreconditionerDeleter<ILU>(ilu);
     }
+#endif
 
 //! \brief Creates and initializes a unique pointer to an sequential ILU0 preconditioner.
 //! \param A     The matrix of the linear system to solve.
@@ -164,6 +166,7 @@ createILUnPtr(const M& A, int ilu_n, double relax, const Dune::Amg::SequentialIn
     return std::shared_ptr<Dune::SeqILUn<M,X,X> >(new Dune::SeqILUn<M,X,X>( A, ilu_n, relax) );
 }
 
+#if HAVE_MPI
 template<class ILU, class I1, class I2>
 struct SelectParallelILUSharedPtr
 {
@@ -218,6 +221,7 @@ createILUnPtr(const M& A, int ilu_n, double 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.
@@ -230,6 +234,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
     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.
@@ -251,6 +256,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
     return EllipticPreconditionerPointer(new ParallelPreconditioner(*ilu, comm),
                                          createParallelDeleter(*ilu, comm));
 }
+#endif
 } // end namespace
 
 

From dd63c2489f4b7f9577c1e264f1ff617df90b69c9 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Fri, 6 Feb 2015 10:41:10 +0100
Subject: [PATCH 09/15] Do not reimplement null_deleter from dune-common.

Instead we use Dune::stackobject_to_shared_ptr to create a
shared_ptr that does not delete the pointer.
---
 .../ExtractParallelGridInformationToISTL.cpp    | 17 +++++------------
 1 file changed, 5 insertions(+), 12 deletions(-)

diff --git a/opm/autodiff/ExtractParallelGridInformationToISTL.cpp b/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
index 42f30a220..7e09d7b50 100644
--- a/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
+++ b/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
@@ -21,29 +21,22 @@
 #include "ExtractParallelGridInformationToISTL.hpp"
 #include <opm/core/linalg/ParallelIstlInformation.hpp>
 #include <dune/common/version.hh>
-
+#include <dune/common/shared_ptr.hh>
  namespace Opm
 {
 #if defined(HAVE_DUNE_CORNERPOINT) && DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
 #if defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL)
 // Extracts the information about the data decomposition from the grid for dune-istl
-struct NullDeleter
-{
-    void operator()(void*)
-    {}
-};
-
 void extractParallelGridInformationToISTL(boost::any& anyComm, const Dune::CpGrid& grid)
 {
     if(grid.comm().size()>1)
     {
         // this is a parallel run with distributed data.
-        NullDeleter null_delete;
         Dune::CpGrid& mgrid=const_cast<Dune::CpGrid&>(grid);
-        Dune::CpGrid::ParallelIndexSet* idx=&(mgrid.getCellIndexSet());
-        Dune::CpGrid::RemoteIndices* ridx=&(mgrid.getCellRemoteIndices());
-        anyComm=boost::any(Opm::ParallelISTLInformation(std::shared_ptr<Dune::CpGrid::ParallelIndexSet>(idx, null_delete),
-                                                        std::shared_ptr<Dune::CpGrid::RemoteIndices>(ridx, null_delete),
+        Dune::CpGrid::ParallelIndexSet& idx=mgrid.getCellIndexSet();
+        Dune::CpGrid::RemoteIndices& ridx=mgrid.getCellRemoteIndices();
+        anyComm=boost::any(Opm::ParallelISTLInformation(Dune::stackobject_to_shared_ptr(idx),
+                                                        Dune::stackobject_to_shared_ptr(ridx),
                                                         grid.comm()));
     }
 }

From bf281a74b8ada4398140f3f254149dfa64f5b12c Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Fri, 6 Feb 2015 10:42:19 +0100
Subject: [PATCH 10/15] Make symbol extractParallelGridInformationToISTL
 available if CpGrid is there.

---
 opm/autodiff/ExtractParallelGridInformationToISTL.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/opm/autodiff/ExtractParallelGridInformationToISTL.cpp b/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
index 7e09d7b50..f5e763cdd 100644
--- a/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
+++ b/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
@@ -24,8 +24,8 @@
 #include <dune/common/shared_ptr.hh>
  namespace Opm
 {
-#if defined(HAVE_DUNE_CORNERPOINT) && DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
-#if defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL)
+#if defined(HAVE_DUNE_CORNERPOINT)
+#if defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL) && DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
 // Extracts the information about the data decomposition from the grid for dune-istl
 void extractParallelGridInformationToISTL(boost::any& anyComm, const Dune::CpGrid& grid)
 {
@@ -44,6 +44,6 @@ void extractParallelGridInformationToISTL(boost::any& anyComm, const Dune::CpGri
 // Missing support for MPI or dune-istl -> do nothing.
 void extractParallelGridInformationToISTL(boost::any&, const Dune::CpGrid&)
 {}
-#endif //defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL)
+#endif //defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL)  && DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
 #endif //defined(HAVE_DUNE_CORNERPOINT)
 } // end namespace Opm

From 31da0de90920a6b72a21d0ed16dcded80538e01e Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Wed, 11 Feb 2015 10:53:37 +0100
Subject: [PATCH 11/15] Fixes indentation in CMakeLists_files.cmake

---
 CMakeLists_files.cmake | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 9ebb3a6c4..4cf2ce34d 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -102,7 +102,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/CPRPreconditioner.hpp
 	opm/autodiff/fastSparseProduct.hpp
 	opm/autodiff/DuneMatrix.hpp
-  opm/autodiff/ExtractParallelGridInformationToISTL.hpp
+	opm/autodiff/ExtractParallelGridInformationToISTL.hpp
 	opm/autodiff/GeoProps.hpp
 	opm/autodiff/GridHelpers.hpp
 	opm/autodiff/ImpesTPFAAD.hpp

From a37421d2adb041705955d1a75caadcee4083cb04 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Wed, 11 Feb 2015 10:58:59 +0100
Subject: [PATCH 12/15] Rely on delete being null safe.

---
 opm/autodiff/CPRPreconditioner.hpp | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp
index 1d078345a..22957e0e5 100644
--- a/opm/autodiff/CPRPreconditioner.hpp
+++ b/opm/autodiff/CPRPreconditioner.hpp
@@ -68,11 +68,8 @@ public:
     template<class T>
     void operator()(T* pt)
     {
-        if(ilu_)
-        {
-            delete ilu_;
-        }
         delete pt;
+        delete ilu_;
     }
 private:
     PREC* ilu_;

From 11848acf0c8627f8e66b5ea4bb5a08f11237d160 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Wed, 11 Feb 2015 11:01:34 +0100
Subject: [PATCH 13/15] Added braces around else statement (coding guidelines)

---
 opm/autodiff/CPRPreconditioner.hpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp
index 22957e0e5..3a58b0013 100644
--- a/opm/autodiff/CPRPreconditioner.hpp
+++ b/opm/autodiff/CPRPreconditioner.hpp
@@ -522,7 +522,9 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
               amg_ = std::unique_ptr< AMG > (new AMG(*opAe_, criterion, smootherArgs));
             }
             else
+            {
                 precond_ = createEllipticPreconditionerPointer<M,X>( Ae_, relax_, comm);
+            }
        }
     };
 

From be7221aa9beda94c52532c89715497698cedb43b Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Wed, 11 Feb 2015 11:22:27 +0100
Subject: [PATCH 14/15] Moved output parameters to the end of the list in
 extractParallelGridInformationToISTL

---
 examples/sim_fibo_ad_cp.cpp                           | 2 +-
 opm/autodiff/ExtractParallelGridInformationToISTL.cpp | 4 ++--
 opm/autodiff/ExtractParallelGridInformationToISTL.hpp | 6 +++---
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/examples/sim_fibo_ad_cp.cpp b/examples/sim_fibo_ad_cp.cpp
index b309727b6..af40e230c 100644
--- a/examples/sim_fibo_ad_cp.cpp
+++ b/examples/sim_fibo_ad_cp.cpp
@@ -213,7 +213,7 @@ try
     std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
     
     boost::any parallel_information;
-    Opm::extractParallelGridInformationToISTL(parallel_information, *grid);
+    Opm::extractParallelGridInformationToISTL(*grid, parallel_information);
     if (param.getDefault("use_cpr", true)) {
         fis_solver.reset(new NewtonIterationBlackoilCPR(param, parallel_information));
     } else {
diff --git a/opm/autodiff/ExtractParallelGridInformationToISTL.cpp b/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
index f5e763cdd..c725bcf8c 100644
--- a/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
+++ b/opm/autodiff/ExtractParallelGridInformationToISTL.cpp
@@ -27,7 +27,7 @@
 #if defined(HAVE_DUNE_CORNERPOINT)
 #if defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL) && DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
 // Extracts the information about the data decomposition from the grid for dune-istl
-void extractParallelGridInformationToISTL(boost::any& anyComm, const Dune::CpGrid& grid)
+void extractParallelGridInformationToISTL(const Dune::CpGrid& grid, boost::any& anyComm)
 {
     if(grid.comm().size()>1)
     {
@@ -42,7 +42,7 @@ void extractParallelGridInformationToISTL(boost::any& anyComm, const Dune::CpGri
 }
 #else
 // Missing support for MPI or dune-istl -> do nothing.
-void extractParallelGridInformationToISTL(boost::any&, const Dune::CpGrid&)
+void extractParallelGridInformationToISTL(const Dune::CpGrid&, boost::any&)
 {}
 #endif //defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL)  && DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
 #endif //defined(HAVE_DUNE_CORNERPOINT)
diff --git a/opm/autodiff/ExtractParallelGridInformationToISTL.hpp b/opm/autodiff/ExtractParallelGridInformationToISTL.hpp
index ab57d456d..71d2a3bb2 100644
--- a/opm/autodiff/ExtractParallelGridInformationToISTL.hpp
+++ b/opm/autodiff/ExtractParallelGridInformationToISTL.hpp
@@ -33,11 +33,11 @@ namespace Opm
 /// about the data decompoisition and convert it to the format expected by the linear algebra
 /// of dune-istl.
 /// \warn if there is no support for dune-istl and MPI then this functio does not do anything.
-/// \param anyComm The handle to to store the information in. If grid is a parallel grid
+/// \param[in] grid The grid to inspect.
+/// \param[out] anyComm The handle to to store the information in. If grid is a parallel grid
 /// then this will ecapsulate an instance of ParallelISTLInformation.
-/// \param grid The grid to inspect.
 
-void extractParallelGridInformationToISTL(boost::any& anyComm, const Dune::CpGrid& grid);
+void extractParallelGridInformationToISTL(const Dune::CpGrid& grid, boost::any& anyComm);
 
 } // end namespace Opm
 #endif //defined(HAVE_DUNE_CORNERPOINT)

From d1239d0a35e4a103dab66dc0bcaf4016fbdc1bf5 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Wed, 11 Feb 2015 11:23:43 +0100
Subject: [PATCH 15/15] Adds space between parameters and stops eliding zeros.

---
 opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index 8a2883699..da192aed6 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -1908,7 +1908,7 @@ namespace {
                 boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
             // Compute the global number of cells and porevolume
             std::vector<int> v(nc, 1);
-            auto nc_and_pv = std::tuple<int, double>(0,.0);
+            auto nc_and_pv = std::tuple<int, double>(0, 0.0);
             auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
                                                         Opm::Reduction::makeGlobalSumFunctor<double>());
             auto nc_and_pv_containers  = std::make_tuple(v, geo_.poreVolume());
@@ -1917,21 +1917,21 @@ namespace {
             for ( int idx=0; idx<MaxNumPhases; ++idx )
             {
                 if (active_[idx]) {
-                    auto values     = std::tuple<double,double,double>(0.,0.,0.);
+                    auto values     = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
                     auto containers = std::make_tuple(B.col(idx),
                                                       tempV.col(idx),
                                                       R.col(idx));
                     auto operators  = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
                                                       Opm::Reduction::makeGlobalMaxFunctor<double>(),
                                                       Opm::Reduction::makeGlobalSumFunctor<double>());
-                    info.computeReduction(containers,operators,values);
+                    info.computeReduction(containers, operators, values);
                     B_avg[idx]    = std::get<0>(values)/std::get<0>(nc_and_pv);
                     maxCoeff[idx] = std::get<1>(values);
                     R_sum[idx]    = std::get<2>(values);
                 }
                 else
                 {
-                    R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.;
+                    R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.0;
                 }
             }
             // Compute pore volume
@@ -1949,7 +1949,7 @@ namespace {
                 }
                 else
                 {
-                    R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.;
+                    R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.0;
                 }
             }
             // Compute total pore volume