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

@@ -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 */

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_;
         };
     }
 }

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