From 5fce54fd447a6e9412c5c9aaebcb5873f6511597 Mon Sep 17 00:00:00 2001
From: Robert Kloefkorn <robert.kloefkorn@iris.no>
Date: Thu, 17 Nov 2016 16:17:41 +0100
Subject: [PATCH] [bugfix][ISTLSolver] make code compile with AMG when
 different matrix operator is used.

---
 opm/autodiff/ISTLSolver.hpp | 37 ++++++++++++++++++++++++++++---------
 1 file changed, 28 insertions(+), 9 deletions(-)

diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp
index b7ca977b1..a5132c4a1 100644
--- a/opm/autodiff/ISTLSolver.hpp
+++ b/opm/autodiff/ISTLSolver.hpp
@@ -191,8 +191,8 @@ namespace Opm
         /// \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 POrComm>
-        void constructPreconditionerAndSolve(O& opA,
+        template<int category=Dune::SolverCategory::sequential, class LinearOperator, class POrComm>
+        void constructPreconditionerAndSolve(LinearOperator& linearOperator,
                                              Vector& x, Vector& istlb,
                                              const POrComm& parallelInformation_arg,
                                              Dune::InverseOperatorResult& result) const
@@ -211,21 +211,33 @@ namespace Opm
             {
                 typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm>  CPRSelectorType;
                 typedef typename CPRSelectorType::AMG AMG;
+                typedef typename CPRSelectorType::Operator MatrixOperator;
+
                 std::unique_ptr< AMG > amg;
+                std::unique_ptr< MatrixOperator > opA;
+
+                if( ! std::is_same< LinearOperator, MatrixOperator > :: value )
+                {
+                    // create new operator in case linear operator and matrix operator differ
+                    opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) );
+                }
+
+                const double relax = 1.0;
+
                 // Construct preconditioner.
-                constructAMGPrecond(opA, parallelInformation_arg, amg);
+                constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
 
                 // Solve.
-                solve(opA, x, istlb, *sp, *amg, result);
+                solve(linearOperator, x, istlb, *sp, *amg, result);
             }
             else
 #endif
             {
                 // Construct preconditioner.
-                auto precond = constructPrecond(opA, parallelInformation_arg);
+                auto precond = constructPrecond(linearOperator, parallelInformation_arg);
 
                 // Solve.
-                solve(opA, x, istlb, *sp, *precond, result);
+                solve(linearOperator, x, istlb, *sp, *precond, result);
             }
         }
 
@@ -252,11 +264,18 @@ namespace Opm
         }
 #endif
 
-        template <class Operator, class POrComm, class AMG >
+        template <class LinearOperator, class MatrixOperator, class POrComm, class AMG >
         void
-        constructAMGPrecond(Operator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg ) const
+        constructAMGPrecond(LinearOperator& linearOperator, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
+        {
+            ISTLUtility::createAMGPreconditionerPointer( *opA, relax, comm, amg );
+        }
+
+
+        template <class MatrixOperator, class POrComm, class AMG >
+        void
+        constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
         {
-            const double relax = 1.0;
             ISTLUtility::createAMGPreconditionerPointer( opA, relax, comm, amg );
         }