diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp
index 1e8df834..cd6b996b 100644
--- a/examples/sim_wateroil.cpp
+++ b/examples/sim_wateroil.cpp
@@ -363,7 +363,7 @@ main(int argc, char** argv)
     int num_wells = 0;
     if (wells->c_wells()) {
         num_wells = wells->c_wells()->number_of_wells;
-        well_state.init(wells->c_wells());
+        well_state.init(wells->c_wells(), state);
         well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(num_wells), 0.0);
         wellreport.push(*props, *wells->c_wells(),
                         state.pressure(), state.surfacevol(), state.saturation(),
diff --git a/opm/core/WellState.hpp b/opm/core/WellState.hpp
index 442b7533..56ede664 100644
--- a/opm/core/WellState.hpp
+++ b/opm/core/WellState.hpp
@@ -29,11 +29,18 @@ namespace Opm
     class WellState
     {
     public:
-        void init(const Wells* wells)
+        template <class State>
+        void init(const Wells* wells, const State& state)
         {
             if (wells) {
-                bhp_.resize(wells->number_of_wells);
-                perfrates_.resize(wells->well_connpos[wells->number_of_wells]);
+                const int nw = wells->number_of_wells;
+                bhp_.resize(nw);
+                // Initialize bhp to be pressure in first perforation cell.
+                for (int w = 0; w < nw; ++w) {
+                    const int cell = wells->well_cells[wells->well_connpos[w]];
+                    bhp_[w] = state.pressure()[cell];
+                }
+                perfrates_.resize(wells->well_connpos[nw]);
             }
         }
 
diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp
index 52acc35f..0ee8bf8b 100644
--- a/opm/core/pressure/CompressibleTpfa.cpp
+++ b/opm/core/pressure/CompressibleTpfa.cpp
@@ -44,7 +44,10 @@ namespace Opm
     /// \param[in] grid          A 2d or 3d grid.
     /// \param[in] props         Rock and fluid properties.
     /// \param[in] linsolver     Linear solver to use.
-    /// \param[in] gravity       Gravity vector. If nonzero, the array should
+    /// \param[in] residual_tol  Solution accepted if inf-norm of residual is smaller.
+    /// \param[in] change_tol    Solution accepted if inf-norm of change is smaller.
+    /// \param[in] maxiter       Maximum acceptable 
+    /// \param[in] gravity       Gravity vector. If non-null, the array should
     ///                          have D elements.
     /// \param[in] wells         The wells argument. Will be used in solution, 
     ///                          is ignored if NULL.
@@ -111,6 +114,7 @@ namespace Opm
                                  WellState& well_state)
     {
         const int nc = grid_.number_of_cells;
+        const int nw = wells_->number_of_wells;
 
         // Set up dynamic data.
         computePerSolveDynamicData(dt, state, well_state);
@@ -124,20 +128,29 @@ namespace Opm
         double res_norm = residualNorm();
         std::cout << "\nIteration         Residual        Change in p\n"
                   << std::setw(9) << iter
-                  << std::setw(18) << res_norm << std::endl;
-        for (; (iter < maxiter_) && (res_norm > residual_tol_); ++iter) {
+                  << std::setw(18) << res_norm
+                  << std::setw(18) << '*' << std::endl;
+        for (; (iter < maxiter_) && (res_norm > residual_tol_); ) {
             // Solve for increment in Newton method:
             //   incr = x_{n+1} - x_{n} = -J^{-1}F
             // (J is Jacobian matrix, F is residual)
             solveIncrement();
+            ++iter;
 
             // Update pressure vars with increment.
-            std::copy(pressure_increment_.begin(), pressure_increment_.begin() + nc, state.pressure().begin());
-            std::copy(pressure_increment_.begin() + nc, pressure_increment_.end(), well_state.bhp().begin());
+            for (int c = 0; c < nc; ++c) {
+                state.pressure()[c] += pressure_increment_[c];
+            }
+            for (int w = 0; w < nw; ++w) {
+                well_state.bhp()[w] += pressure_increment_[nc + w];
+            }
 
             // Stop iterating if increment is small.
             inc_norm = incrementNorm();
             if (inc_norm <= change_tol_) {
+                std::cout << std::setw(9) << iter
+                          << std::setw(18) << '*' 
+                          << std::setw(18) << inc_norm << std::endl;
                 break;
             }
 
@@ -150,7 +163,7 @@ namespace Opm
             // Update residual norm.
             res_norm = residualNorm();
 
-            std::cout << std::setw(9) << iter + 1
+            std::cout << std::setw(9) << iter
                       << std::setw(18) << res_norm 
                       << std::setw(18) << inc_norm << std::endl;
         }
@@ -408,11 +421,9 @@ namespace Opm
         wellperf_phasemob_.resize(nperf*np);
         // The A matrix is set equal to the perforation grid cells'
         // matrix, for both injectors and producers.
-        // The mobilities are all set equal to the total mobility for the
-        // cell for injectors, and equal to individual phase mobilities
-        // for producers.
+        // The mobilities are set equal to the perforation grid cells'
+        // mobilities, for both injectors and producers.
         for (int w = 0; w < nw; ++w) {
-            bool is_injector = wells_->type[w] == INJECTOR;
             for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
                 const int c = wells_->well_cells[j];
                 const double* cA = &cell_A_[np*np*c];
@@ -420,15 +431,7 @@ namespace Opm
                 std::copy(cA, cA + np*np, wpA);
                 const double* cM = &cell_phasemob_[np*c];
                 double* wpM = &wellperf_phasemob_[np*j];
-                if (is_injector) {
-                    double totmob = 0.0;
-                    for (int phase = 0; phase < np; ++phase) {
-                        totmob += cM[phase];
-                    }
-                    std::fill(wpM, wpM + np, totmob);
-                } else {
-                    std::copy(cM, cM + np, wpM);
-                }
+                std::copy(cM, cM + np, wpM);
             }
         }
     }