From a011244d9caf10d5422ecd6d5a63269f071e2a95 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Thu, 23 Jul 2020 10:31:37 +0200
Subject: [PATCH] refactored invDxDirect to allow for reusing the
 decomposition.

---
 opm/simulators/wells/MSWellHelpers.hpp | 23 +++++++++++++++++------
 1 file changed, 17 insertions(+), 6 deletions(-)

diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp
index a3b0cfcc2..f92d6b732 100644
--- a/opm/simulators/wells/MSWellHelpers.hpp
+++ b/opm/simulators/wells/MSWellHelpers.hpp
@@ -1,6 +1,7 @@
 /*
   Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
   Copyright 2017 Statoil ASA.
+  Copyright 2020 Equinor ASA.
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -36,17 +37,17 @@ namespace Opm {
 
 namespace mswellhelpers
 {
-    // obtain y = D^-1 * x with a direct solver
+
+#if HAVE_UMFPACK
+    /// Applies umfpack and checks for singularity
     template <typename MatrixType, typename VectorType>
     VectorType
-    invDXDirect(const MatrixType& D, VectorType x)
+    applyUMFPack(Dune::UMFPack<MatrixType>& linsolver, VectorType x)
     {
-#if HAVE_UMFPACK
+        // The copy of x seems mandatory for calling UMFPack!
         VectorType y(x.size());
         y = 0.;
 
-        Dune::UMFPack<MatrixType> linsolver(D, 0);
-
         // Object storing some statistics about the solving process
         Dune::InverseOperatorResult res;
 
@@ -62,8 +63,18 @@ namespace mswellhelpers
                 }
             }
         }
-
         return y;
+    }
+#endif
+
+    // obtain y = D^-1 * x with a direct solver
+    template <typename MatrixType, typename VectorType>
+    VectorType
+    invDXDirect(const MatrixType& D, VectorType x)
+    {
+#if HAVE_UMFPACK
+        Dune::UMFPack<MatrixType> linsolver(D, 0);
+        return applyUMFPack(linsolver, x);
 #else
         // this is not thread safe
         OPM_THROW(std::runtime_error, "Cannot use invDXDirect() without UMFPACK. "