From d5010efdba1f276c6f82e656f76d33f831adadaf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= <Bard.Skaflestad@sintef.no>
Date: Mon, 3 Oct 2011 11:21:41 +0200
Subject: [PATCH] Solvers need to be fully formed objects.  Store references.

For instance, an iterative linear solver is likely to need some kind
of stopping criterion and the implicit solver should not be tasked
with defining the criterion.
---
 src/ImplicitTransport.hpp         |  7 ++++++-
 src/JacobianSystem.hpp            | 21 ++++++++++++---------
 src/SinglePointUpwindTwoPhase.hpp |  7 +++++--
 3 files changed, 23 insertions(+), 12 deletions(-)

diff --git a/src/ImplicitTransport.hpp b/src/ImplicitTransport.hpp
index feca16fc..d68270de 100644
--- a/src/ImplicitTransport.hpp
+++ b/src/ImplicitTransport.hpp
@@ -59,6 +59,8 @@ namespace Opm {
               class JacobianSystem>
     class ImplicitTransport : private Model {
     public:
+        ImplicitTransport(JacobianSystem& sys) : sys_(sys) {}
+
         template <class Grid          ,
                   class SourceTerms   ,
                   class ReservoirState>
@@ -112,13 +114,16 @@ namespace Opm {
         }
 
     private:
+        ImplicitTransport           (const ImplicitTransport&);
+        ImplicitTransport& operator=(const ImplicitTransport&);
+
         using Model::initStep;
         using Model::initIteration;
         using Model::finishIteration;
         using Model::finishStep;
 
         ImplicitAssembly<Model> asm_;
-        JacobianSystem          sys_;
+        JacobianSystem&         sys_;
     };
 }
 #endif  /* OPM_IMPLICITTRANSPORT_HPP_HEADER */
diff --git a/src/JacobianSystem.hpp b/src/JacobianSystem.hpp
index c19b7f77..b7c9cb7f 100644
--- a/src/JacobianSystem.hpp
+++ b/src/JacobianSystem.hpp
@@ -203,12 +203,11 @@ namespace Opm {
             matrix();
         } */;
 
-        template <class Matrix, class Vector,
-                  template <class Mat,
-                            class Vec>
-                  class LinSolve>
+        template <class Matrix, class Vector, class LinSolve>
         class JacobianSystem {
         public:
+            JacobianSystem(LinSolve& solver) : solver_(solver) {}
+
             typedef Matrix matrix_type;
 
             MatrixBlockAssembler<Matrix>& matrix() { return mba_ ; }
@@ -216,14 +215,18 @@ namespace Opm {
 
             void
             solve() {
-                LinSolve<Matrix,Vector> ls(mba_.matrix());
-
-                ls.solve(vecs_.residual(), vecs_.increment());
+                solver_.solve(mba_ .matrix   (),
+                              vecs_.residual (),
+                              vecs_.increment());
             }
 
         private:
-            MatrixBlockAssembler<Matrix> mba_ ;
-            NewtonVectors       <Vector> vecs_;
+            JacobianSystem           (const JacobianSystem&);
+            JacobianSystem& operator=(const JacobianSystem&);
+
+            MatrixBlockAssembler<Matrix> mba_   ;
+            NewtonVectors       <Vector> vecs_  ;
+            LinSolve&                    solver_;
         };
     }
 }
diff --git a/src/SinglePointUpwindTwoPhase.hpp b/src/SinglePointUpwindTwoPhase.hpp
index 32608802..b444a275 100644
--- a/src/SinglePointUpwindTwoPhase.hpp
+++ b/src/SinglePointUpwindTwoPhase.hpp
@@ -258,11 +258,14 @@ namespace Opm {
         initIteration(const ReservoirState& state,
                       const Grid&           g    ,
                       JacobianSystem&       sys  ) {
+
+            std::array<double, 2  > s   ;
             std::array<double, 2*2> dmob;
 
-            const double *s = &state.saturation[0*2 + 0];
+            for (int c = 0; c < g.number_of_cells; ++c) {
+                s[0] = state.saturation[c*2 + 0] + sys.vector().solution()[c];
+                s[1] = 1 - s[0];
 
-            for (int c = 0; c < g.number_of_cells; ++c, s += 0) {
                 this->mobility(c, s, store_.mob(c), dmob);
 
                 store_.dmob(c)[0] =  dmob[0*2 + 0];