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];