From 8e78f022303efd08e2b493c257b19474ec77e28e Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Mon, 26 Feb 2018 12:16:04 +0100
Subject: [PATCH] Use Dune::FMatrixHelp::invertMatrix for inversion

Thus we always take advantage of the specializations in ISTLSolver.hpp
and do need to take care about the type of the matrix block.
---
 opm/autodiff/StandardWell.hpp      | 8 +-------
 opm/autodiff/StandardWell_impl.hpp | 5 ++++-
 2 files changed, 5 insertions(+), 8 deletions(-)

diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp
index 8378b6b63..ac1e18502 100644
--- a/opm/autodiff/StandardWell.hpp
+++ b/opm/autodiff/StandardWell.hpp
@@ -27,6 +27,7 @@
 #include <opm/autodiff/WellInterface.hpp>
 #include <opm/autodiff/ISTLSolver.hpp>
 #include <opm/autodiff/RateConverter.hpp>
+#include <opm/autodiff/ISTLSolver.hpp>
 
 namespace Opm
 {
@@ -88,15 +89,8 @@ namespace Opm
         typedef Dune::FieldVector<Scalar, numWellEq> VectorBlockWellType;
         typedef Dune::BlockVector<VectorBlockWellType> BVectorWell;
 
-#if  DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2 , 5, 1)
-        // 3x3 matrix block inversion was unstable from at least 2.3 until and
-        // including 2.5.0
         // the matrix type for the diagonal matrix D
         typedef Dune::FieldMatrix<Scalar, numWellEq, numWellEq > DiagMatrixBlockWellType;
-#else
-        // the matrix type for the diagonal matrix D
-        typedef Dune::MatrixBlock<Scalar, numWellEq, numWellEq > DiagMatrixBlockWellType;
-#endif
 
         typedef Dune::BCRSMatrix <DiagMatrixBlockWellType> DiagMatWell;
 
diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp
index 18b71e3f1..8d82980d8 100644
--- a/opm/autodiff/StandardWell_impl.hpp
+++ b/opm/autodiff/StandardWell_impl.hpp
@@ -688,7 +688,10 @@ namespace Detail
         }
 
         // do the local inversion of D.
-        invDuneD_[0][0].invert();
+        // we do this manually with invertMatrix to always get our
+        // specializations in for 3x3 and 4x4 matrices.
+        auto original = invDuneD_[0][0];
+        Dune::FMatrixHelp::invertMatrix(original, invDuneD_[0][0]);
 
         if ( param_.matrix_add_well_contributions_ )
         {