diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 523f4a92f..44698d4ed 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -35,8 +35,8 @@ list (APPEND MAIN_SOURCE_FILES
 	opm/autodiff/GridHelpers.cpp
 	opm/autodiff/ImpesTPFAAD.cpp
 	opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
-#	opm/autodiff/SimulatorIncompTwophaseAd.cpp
-#	opm/autodiff/TransportSolverTwophaseAd.cpp
+	opm/autodiff/SimulatorIncompTwophaseAd.cpp
+	opm/autodiff/TransportSolverTwophaseAd.cpp
 	opm/autodiff/BlackoilPropsAdFromDeck.cpp
 #	opm/autodiff/SolventPropsAdFromDeck.cpp
 	opm/autodiff/BlackoilModelParameters.cpp
diff --git a/opm/autodiff/AutoDiffMatrix.hpp b/opm/autodiff/AutoDiffMatrix.hpp
index 646b48e4a..b3b4131d2 100644
--- a/opm/autodiff/AutoDiffMatrix.hpp
+++ b/opm/autodiff/AutoDiffMatrix.hpp
@@ -265,6 +265,40 @@ namespace Opm
 
 
 
+        AutoDiffMatrix operator/(const double rhs) const
+        {
+            switch (type_) {
+            case Z:
+                return *this;
+            case I:
+                {
+                    AutoDiffMatrix retval(*this);
+                    retval.type_ = D;
+                    retval.d_.assign(rows_, 1.0/rhs);
+                    return retval;
+                }
+            case D:
+                {
+                    AutoDiffMatrix retval(*this);
+                    for (double& elem : retval.d_) {
+                        elem /= rhs;
+                    }
+                    return retval;
+                }
+            case S:
+                {
+                    AutoDiffMatrix retval(*this);
+                    retval.s_ /= rhs;
+                    return retval;
+                }
+            }
+        }
+
+
+
+
+
+
         Eigen::VectorXd operator*(const Eigen::VectorXd& rhs) const
         {
             assert(cols_ == rhs.size());
diff --git a/opm/autodiff/TransportSolverTwophaseAd.cpp b/opm/autodiff/TransportSolverTwophaseAd.cpp
index a6ab83f71..4d457a848 100644
--- a/opm/autodiff/TransportSolverTwophaseAd.cpp
+++ b/opm/autodiff/TransportSolverTwophaseAd.cpp
@@ -96,6 +96,7 @@ namespace Opm
         {
             typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
             typedef Eigen::Array<double, Eigen::Dynamic, 4, Eigen::RowMajor> FourCol;
+            typedef Eigen::SparseMatrix<double> S;
             typedef typename ADB::V V;
             typedef typename ADB::M M;
             const int nc = props.numCells();
@@ -117,8 +118,8 @@ namespace Opm
             //    dkro/dsw = col(1) - col(3)
             V dkrw = dkr.leftCols<1>() - dkr.rightCols<2>().leftCols<1>();
             V dkro = dkr.leftCols<2>().rightCols<1>() - dkr.rightCols<1>();
-            M krwjac(nc,nc);
-            M krojac(nc,nc);
+            S krwjac(nc,nc);
+            S krojac(nc,nc);
             auto sizes = Eigen::ArrayXi::Ones(nc);
             krwjac.reserve(sizes);
             krojac.reserve(sizes);
@@ -127,8 +128,8 @@ namespace Opm
                 krojac.insert(c,c) = dkro(c);
             }
             const double* mu = props.viscosity();
-            std::vector<M> dmw = { krwjac/mu[0] };
-            std::vector<M> dmo = { krojac/mu[1] };
+            std::vector<M> dmw = { M(krwjac)/mu[0] };
+            std::vector<M> dmo = { M(krojac)/mu[1] };
 
             std::vector<ADB> pmobc = { ADB::function(krw / mu[0], std::move(dmw)) ,
                                        ADB::function(kro / mu[1], std::move(dmo)) };
@@ -229,7 +230,8 @@ namespace Opm
             std::cout << "Residual l2-norm = " << res_norm << std::endl;
 
             // Solve linear system.
-            Eigen::SparseMatrix<double, Eigen::RowMajor> smatr = transport_residual.derivative()[0];
+            Eigen::SparseMatrix<double, Eigen::RowMajor> smatr;
+            transport_residual.derivative()[0].toSparse(smatr);
             assert(smatr.isCompressed());
             V ds(nc);
             LinearSolverInterface::LinearSolverReport rep