incompressible two phase fully implicit simulator.

opm/polymer/fullyimplicit/AutoDiff.hpp

@@ -0,0 +1,302 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2013 Statoil ASA.
+
+  This file is part of the Open Porous Media Project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_AUTODIFF_HPP_HEADER
+#define OPM_AUTODIFF_HPP_HEADER
+
+#include <cmath>
+
+namespace Opm
+{
+
+    /// A simple class for forward-mode automatic differentiation.
+    ///
+    /// The class represents a single value and a single derivative.
+    /// Only basic arithmetic operators and a few functions are
+    /// implemented for it, it is mostly intended for simple
+    /// experimentation.
+    template <typename Scalar>
+    class AutoDiff
+    {
+    public:
+        /// Create an AutoDiff object representing a constant, that
+        /// is, its derivative is zero.
+        static AutoDiff
+        constant(const Scalar x)
+        {
+            return function(x, Scalar(0));
+        }
+
+        /// Create an AutoDiff object representing a primary variable,
+        /// that is, its derivative is one.
+        static AutoDiff
+        variable(const Scalar x)
+        {
+            return function(x, Scalar(1));
+        }
+
+        /// Create an AutoDiff object representing a function value
+        /// and its derivative.
+        static AutoDiff
+        function(const Scalar x, const Scalar dx)
+        {
+            return AutoDiff(x, dx);
+        }
+
+        void
+        operator +=(const Scalar& rhs)
+        {
+            val_ += rhs;
+        }
+
+        void
+        operator +=(const AutoDiff& rhs)
+        {
+            val_ += rhs.val_;
+            der_ += rhs.der_;
+        }
+
+        void
+        operator -=(const Scalar& rhs)
+        {
+            val_ -= rhs;
+        }
+
+        void
+        operator -=(const AutoDiff& rhs)
+        {
+            val_ -= rhs.val_;
+            der_ -= rhs.der_;
+        }
+
+        void
+        operator *=(const Scalar& rhs)
+        {
+            val_ *= rhs;
+            der_ *= rhs;
+        }
+
+        void
+        operator *=(const AutoDiff& rhs)
+        {
+            der_   = der_*rhs.val_ + val_*rhs.der_;
+            val_  *= rhs.val_;
+        }
+
+        void
+        operator /=(const Scalar& rhs)
+        {
+            val_ /= rhs;
+            der_ /= rhs;
+        }
+
+        void
+        operator /=(const AutoDiff& rhs)
+        {
+            der_   = (der_*rhs.val_ - val_*rhs.der_) / (rhs.val_ * rhs.val_);
+            val_  /= rhs.val_;
+        }
+
+        template <class Ostream>
+        Ostream&
+        print(Ostream& os) const
+        {
+            os << "(x,dx) = (" << val_ << ',' << der_ << ")";
+
+            return os;
+        }
+
+        const Scalar val() const { return val_; }
+        const Scalar der() const { return der_; }
+
+    private:
+        AutoDiff(const Scalar x, const Scalar dx)
+            : val_(x), der_(dx)
+        {}
+
+        Scalar val_;
+        Scalar der_;
+    };
+
+
+    template <class Ostream, typename Scalar>
+    Ostream&
+    operator<<(Ostream& os, const AutoDiff<Scalar>& fw)
+    {
+        return fw.print(os);
+    }
+
+    template <typename Scalar>
+    AutoDiff<Scalar>
+    operator +(const AutoDiff<Scalar>& lhs,
+               const AutoDiff<Scalar>& rhs)
+    {
+        AutoDiff<Scalar> ret = lhs;
+        ret += rhs;
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    AutoDiff<Scalar>
+    operator +(const T                lhs,
+               const AutoDiff<Scalar>& rhs)
+    {
+        AutoDiff<Scalar> ret = rhs;
+        ret += Scalar(lhs);
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    AutoDiff<Scalar>
+    operator +(const AutoDiff<Scalar>& lhs,
+               const T                rhs)
+    {
+        AutoDiff<Scalar> ret = lhs;
+        ret += Scalar(rhs);
+
+        return ret;
+    }
+
+    template <typename Scalar>
+    AutoDiff<Scalar>
+    operator -(const AutoDiff<Scalar>& lhs,
+               const AutoDiff<Scalar>& rhs)
+    {
+        AutoDiff<Scalar> ret = lhs;
+        ret -= rhs;
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    AutoDiff<Scalar>
+    operator -(const T                lhs,
+               const AutoDiff<Scalar>& rhs)
+    {
+        return AutoDiff<Scalar>::function(Scalar(lhs) - rhs.val(), -rhs.der());
+    }
+
+    template <typename Scalar, typename T>
+    AutoDiff<Scalar>
+    operator -(const AutoDiff<Scalar>& lhs,
+               const T                rhs)
+    {
+        AutoDiff<Scalar> ret = lhs;
+        ret -= Scalar(rhs);
+
+        return ret;
+    }
+
+    template <typename Scalar>
+    AutoDiff<Scalar>
+    operator *(const AutoDiff<Scalar>& lhs,
+               const AutoDiff<Scalar>& rhs)
+    {
+        AutoDiff<Scalar> ret = lhs;
+        ret *= rhs;
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    AutoDiff<Scalar>
+    operator *(const T                lhs,
+               const AutoDiff<Scalar>& rhs)
+    {
+        AutoDiff<Scalar> ret = rhs;
+        ret *= Scalar(lhs);
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    AutoDiff<Scalar>
+    operator *(const AutoDiff<Scalar>& lhs,
+               const T                rhs)
+    {
+        AutoDiff<Scalar> ret = lhs;
+        ret *= Scalar(rhs);
+
+        return ret;
+    }
+
+    template <typename Scalar>
+    AutoDiff<Scalar>
+    operator /(const AutoDiff<Scalar>& lhs,
+               const AutoDiff<Scalar>& rhs)
+    {
+        AutoDiff<Scalar> ret = lhs;
+        ret /= rhs;
+
+        return ret;
+    }
+
+    template <typename Scalar, typename T>
+    AutoDiff<Scalar>
+    operator /(const T                lhs,
+               const AutoDiff<Scalar>& rhs)
+    {
+        Scalar a =  Scalar(lhs) / rhs.val();
+        Scalar b = -Scalar(lhs) / (rhs.val() * rhs.val());
+
+        return AutoDiff<Scalar>::function(a, b);
+    }
+
+    template <typename Scalar, typename T>
+    AutoDiff<Scalar>
+    operator /(const AutoDiff<Scalar>& lhs,
+               const T                rhs)
+    {
+        Scalar a = lhs.val() / Scalar(rhs);
+        Scalar b = lhs.der() / Scalar(rhs);
+
+        return AutoDiff<Scalar>::function(a, b);
+    }
+
+    template <typename Scalar>
+    AutoDiff<Scalar>
+    cos(const AutoDiff<Scalar>& x)
+    {
+        Scalar a = std::cos(x.val());
+        Scalar b = -std::sin(x.val()) * x.der();
+
+        return AutoDiff<Scalar>::function(a, b);
+    }
+
+    template <typename Scalar>
+    AutoDiff<Scalar>
+    sqrt(const AutoDiff<Scalar>& x)
+    {
+        Scalar a = std::sqrt(x.val());
+        Scalar b = (Scalar(1.0) / (Scalar(2.0) * a)) * x.der();
+
+        return AutoDiff<Scalar>::function(a, b);
+    }
+
+} // namespace Opm
+
+namespace std {
+    using Opm::cos;
+    using Opm::sqrt;
+}
+
+#endif  /* OPM_AUTODIFF_HPP_HEADER */

opm/polymer/fullyimplicit/AutoDiffBlock.hpp

@@ -0,0 +1,506 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_AUTODIFFBLOCK_HEADER_INCLUDED
+#define OPM_AUTODIFFBLOCK_HEADER_INCLUDED
+
+#include <Eigen/Eigen>
+#include <Eigen/Sparse>
+#include <vector>
+#include <cassert>
+
+namespace Opm
+{
+
+    /// A class for forward-mode automatic differentiation with vector
+    /// values and sparse jacobian matrices.
+    ///
+    /// The class contains a (column) vector of values and multiple
+    /// sparse matrices representing its partial derivatives. Each
+    /// such matrix has a number of rows equal to the number of rows
+    /// in the value vector, and a number of columns equal to the
+    /// number of discrete variables we want to compute the
+    /// derivatives with respect to. The reason to have multiple such
+    /// jacobians instead of just one is to allow simpler grouping of
+    /// variables, making it easier to implement various
+    /// preconditioning schemes. Only basic arithmetic operators are
+    /// implemented for this class, reflecting our needs so far.
+    ///
+    /// The class is built on the Eigen library, using an Eigen array
+    /// type to contain the values and Eigen sparse matrices for the
+    /// jacobians. The overloaded operators are intended to behave in
+    /// a similar way to Eigen arrays, meaning for example that the *
+    /// operator is elementwise multiplication. The only exception is
+    /// multiplication with a sparse matrix from the left, which is
+    /// treated as an Eigen matrix operation.
+    ///
+    /// There are no public constructors, instead we use the Named
+    /// Constructor pattern. In general, one needs to know which
+    /// variables one wants to compute the derivatives with respect to
+    /// before constructing an AutoDiffBlock. Some of the constructors
+    /// require you to pass a block pattern. This should be a vector
+    /// containing the number of columns you want for each jacobian
+    /// matrix.
+    ///
+    /// For example: you want the derivatives with respect to three
+    /// different variables p, r and s. Assuming that there are 10
+    /// elements in p, and 20 in each of r and s, the block pattern is
+    /// { 10, 20, 20 }. When creating the variables p, r and s in your
+    /// program you have two options:
+    ///     - Use the variable() constructor three times, passing the
+    ///       index (0 for p, 1 for r and 2 for s), initial value of
+    ///       each variable and the block pattern.
+    ///     - Use the variables() constructor passing only the initial
+    ///       values of each variable. The block pattern will be
+    ///       inferred from the size of the initial value vectors.
+    ///       This is usually the simplest option if you have multiple
+    ///       variables. Note that this constructor returns a vector
+    ///       of all three variables, so you need to use index access
+    ///       (operator[]) to get the individual variables (that is p,
+    ///       r and s).
+    ///
+    /// After this, the r variable for example will have a size() of
+    /// 20 and three jacobian matrices. The first is a 20 by 10 zero
+    /// matrix, the second is a 20 by 20 identity matrix, and the
+    /// third is a 20 by 20 zero matrix.
+    template <typename Scalar>
+    class AutoDiffBlock
+    {
+    public:
+        /// Underlying type for values.
+        typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> V;
+        /// Underlying type for jacobians.
+        typedef Eigen::SparseMatrix<Scalar> M;
+
+        /// Construct an empty AutoDiffBlock.
+        static AutoDiffBlock null()
+        {
+            V val;
+            std::vector<M> jac;
+            return AutoDiffBlock(val, jac);
+        }
+
+        /// Create an AutoDiffBlock representing a constant.
+        /// \param[in] val         values 
+        static AutoDiffBlock constant(const V& val)
+        {
+            return AutoDiffBlock(val);
+        }
+
+        /// Create an AutoDiffBlock representing a constant.
+        /// This variant requires specifying the block sizes used
+        /// for the Jacobians even though the Jacobian matrices
+        /// themselves will be zero.
+        /// \param[in] val         values 
+        /// \param[in] blocksizes  block pattern
+        static AutoDiffBlock constant(const V& val, const std::vector<int>& blocksizes)
+        {
+            std::vector<M> jac;
+            const int num_elem = val.size();
+            const int num_blocks = blocksizes.size();
+            // For constants, all jacobian blocks are zero.
+            jac.resize(num_blocks);
+            for (int i = 0; i < num_blocks; ++i) {
+                jac[i] = M(num_elem, blocksizes[i]);
+            }
+            return AutoDiffBlock(val, jac);
+        }
+
+        /// Create an AutoDiffBlock representing a single variable block.
+        /// \param[in] index       index of the variable you are constructing
+        /// \param[in] val         values 
+        /// \param[in] blocksizes  block pattern
+        /// The resulting object will have size() equal to block_pattern[index].
+        /// Its jacobians will all be zero, except for derivative()[index], which
+        /// will be an identity matrix.
+        static AutoDiffBlock variable(const int index, const V& val, const std::vector<int>& blocksizes)
+        {
+            std::vector<M> jac;
+            const int num_elem = val.size();
+            const int num_blocks = blocksizes.size();
+            // First, set all jacobian blocks to zero...
+            jac.resize(num_blocks);
+            for (int i = 0; i < num_blocks; ++i) {
+                jac[i] = M(num_elem, blocksizes[i]);
+            }
+            // ... then set the one corrresponding to this variable to identity.
+            assert(blocksizes[index] == num_elem);
+            jac[index].reserve(Eigen::VectorXi::Constant(val.size(), 1));
+            for (typename M::Index row = 0; row < val.size(); ++row) {
+                jac[index].insert(row, row) = Scalar(1.0);
+            }
+            return AutoDiffBlock(val, jac);
+        }
+
+        /// Create an AutoDiffBlock by directly specifying values and jacobians.
+        /// \param[in] val         values 
+        /// \param[in] jac         vector of jacobians
+        static AutoDiffBlock function(const V& val, const std::vector<M>& jac)
+        {
+            return AutoDiffBlock(val, jac);
+        }
+
+        /// Construct a set of primary variables, each initialized to
+        /// a given vector.
+        static std::vector<AutoDiffBlock> variables(const std::vector<V>& initial_values)
+        {
+            const int num_vars = initial_values.size();
+            std::vector<int> bpat;
+            for (int v = 0; v < num_vars; ++v) {
+                bpat.push_back(initial_values[v].size());
+            }
+            std::vector<AutoDiffBlock> vars;
+            for (int v = 0; v < num_vars; ++v) {
+                vars.emplace_back(variable(v, initial_values[v], bpat));
+            }
+            return vars;
+        }
+
+        /// Elementwise operator +=
+        AutoDiffBlock& operator+=(const AutoDiffBlock& rhs)
+        {
+            if (jac_.empty()) {
+                jac_ = rhs.jac_;
+            } else if (!rhs.jac_.empty()) {
+                assert (numBlocks()    == rhs.numBlocks());
+                assert (value().size() == rhs.value().size());
+
+                const int num_blocks = numBlocks();
+                for (int block = 0; block < num_blocks; ++block) {
+                    assert(jac_[block].rows() == rhs.jac_[block].rows());
+                    assert(jac_[block].cols() == rhs.jac_[block].cols());
+                    jac_[block] += rhs.jac_[block];
+                }
+            }
+
+            val_ += rhs.val_;
+
+            return *this;
+        }
+
+        /// Elementwise operator +
+        AutoDiffBlock operator+(const AutoDiffBlock& rhs) const
+        {
+            if (jac_.empty() && rhs.jac_.empty()) {
+                return constant(val_ + rhs.val_);
+            }
+            if (jac_.empty()) {
+                return val_ + rhs;
+            }
+            if (rhs.jac_.empty()) {
+                return *this + rhs.val_;
+            }
+            std::vector<M> jac = jac_;
+            assert(numBlocks() == rhs.numBlocks());
+            int num_blocks = numBlocks();
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(jac[block].rows() == rhs.jac_[block].rows());
+                assert(jac[block].cols() == rhs.jac_[block].cols());
+                jac[block] += rhs.jac_[block];
+            }
+            return function(val_ + rhs.val_, jac);
+        }
+
+        /// Elementwise operator -
+        AutoDiffBlock operator-(const AutoDiffBlock& rhs) const
+        {
+            if (jac_.empty() && rhs.jac_.empty()) {
+                return constant(val_ - rhs.val_);
+            }
+            if (jac_.empty()) {
+                return val_ - rhs;
+            }
+            if (rhs.jac_.empty()) {
+                return *this - rhs.val_;
+            }
+            std::vector<M> jac = jac_;
+            assert(numBlocks() == rhs.numBlocks());
+            int num_blocks = numBlocks();
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(jac[block].rows() == rhs.jac_[block].rows());
+                assert(jac[block].cols() == rhs.jac_[block].cols());
+                jac[block] -= rhs.jac_[block];
+            }
+            return function(val_ - rhs.val_, jac);
+        }
+
+        /// Elementwise operator *
+        AutoDiffBlock operator*(const AutoDiffBlock& rhs) const
+        {
+            if (jac_.empty() && rhs.jac_.empty()) {
+                return constant(val_ * rhs.val_);
+            }
+            if (jac_.empty()) {
+                return val_ * rhs;
+            }
+            if (rhs.jac_.empty()) {
+                return *this * rhs.val_;
+            }
+            int num_blocks = numBlocks();
+            std::vector<M> jac(num_blocks);
+            assert(numBlocks() == rhs.numBlocks());
+            typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
+            D D1 = val_.matrix().asDiagonal();
+            D D2 = rhs.val_.matrix().asDiagonal();
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(jac_[block].rows() == rhs.jac_[block].rows());
+                assert(jac_[block].cols() == rhs.jac_[block].cols());
+                jac[block] = D2*jac_[block] + D1*rhs.jac_[block];
+            }
+            return function(val_ * rhs.val_, jac);
+        }
+
+        /// Elementwise operator /
+        AutoDiffBlock operator/(const AutoDiffBlock& rhs) const
+        {
+            if (jac_.empty() && rhs.jac_.empty()) {
+                return constant(val_ / rhs.val_);
+            }
+            if (jac_.empty()) {
+                return val_ / rhs;
+            }
+            if (rhs.jac_.empty()) {
+                return *this / rhs.val_;
+            }
+            int num_blocks = numBlocks();
+            std::vector<M> jac(num_blocks);
+            assert(numBlocks() == rhs.numBlocks());
+            typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
+            D D1 = val_.matrix().asDiagonal();
+            D D2 = rhs.val_.matrix().asDiagonal();
+            D D3 = (1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal();
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(jac_[block].rows() == rhs.jac_[block].rows());
+                assert(jac_[block].cols() == rhs.jac_[block].cols());
+                jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]);
+            }
+            return function(val_ / rhs.val_, jac);
+        }
+
+        /// I/O.
+        template <class Ostream>
+        Ostream&
+        print(Ostream& os) const
+        {
+            int num_blocks = jac_.size();
+            os << "Value =\n" << val_ << "\n\nJacobian =\n";
+            for (int i = 0; i < num_blocks; ++i) {
+                os << "Sub Jacobian #" << i << '\n' << jac_[i] << "\n";
+            }
+            return os;
+        }
+
+        /// Number of elements
+        int size() const
+        {
+            return val_.size();
+        }
+
+        /// Number of Jacobian blocks.
+        int numBlocks() const
+        {
+            return jac_.size();
+        }
+
+        /// Sizes (number of columns) of Jacobian blocks.
+        std::vector<int> blockPattern() const
+        {
+            const int nb = numBlocks();
+            std::vector<int> bp(nb);
+            for (int block = 0; block < nb; ++block) {
+                bp[block] = jac_[block].cols();
+            }
+            return bp;
+        }
+
+        /// Function value.
+        const V& value() const
+        {
+            return val_;
+        }
+
+        /// Function derivatives.
+        const std::vector<M>& derivative() const
+        {
+            return jac_;
+        }
+
+    private:
+        AutoDiffBlock(const V& val)
+        : val_(val)
+        {
+        }
+
+        AutoDiffBlock(const V& val,
+                      const std::vector<M>& jac)
+            : val_(val), jac_(jac)
+        {
+#ifndef NDEBUG
+            const int num_elem = val_.size();
+            const int num_blocks = jac_.size();
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(num_elem == jac_[block].rows());
+            }
+#endif
+        }
+
+        V val_;
+        std::vector<M> jac_;
+    };
+
+
+    // ---------  Free functions and operators for AutoDiffBlock  ---------
+
+    /// Stream output.
+    template <class Ostream, typename Scalar>
+    Ostream&
+    operator<<(Ostream& os, const AutoDiffBlock<Scalar>& fw)
+    {
+        return fw.print(os);
+    }
+
+
+    /// Multiply with sparse matrix from the left.
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator*(const typename AutoDiffBlock<Scalar>::M& lhs,
+                                    const AutoDiffBlock<Scalar>& rhs)
+    {
+        int num_blocks = rhs.numBlocks();
+        std::vector<typename AutoDiffBlock<Scalar>::M> jac(num_blocks);
+        assert(lhs.cols() == rhs.value().rows());
+        for (int block = 0; block < num_blocks; ++block) {
+            jac[block] = lhs*rhs.derivative()[block];
+        }
+        typename AutoDiffBlock<Scalar>::V val = lhs*rhs.value().matrix();
+        return AutoDiffBlock<Scalar>::function(val, jac);
+    }
+
+
+    /// Elementwise multiplication with constant on the left.
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator*(const typename AutoDiffBlock<Scalar>::V& lhs,
+                                    const AutoDiffBlock<Scalar>& rhs)
+    {
+        return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) * rhs;
+    }
+
+
+    /// Elementwise multiplication with constant on the right.
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator*(const AutoDiffBlock<Scalar>& lhs,
+                                    const typename AutoDiffBlock<Scalar>::V& rhs)
+    {
+        return rhs * lhs; // Commutative operation.
+    }
+
+
+    /// Elementwise addition with constant on the left.
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator+(const typename AutoDiffBlock<Scalar>::V& lhs,
+                                    const AutoDiffBlock<Scalar>& rhs)
+    {
+        return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) + rhs;
+    }
+
+
+    /// Elementwise addition with constant on the right.
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator+(const AutoDiffBlock<Scalar>& lhs,
+                                    const typename AutoDiffBlock<Scalar>::V& rhs)
+    {
+        return rhs + lhs; // Commutative operation.
+    }
+
+
+    /// Elementwise subtraction with constant on the left.
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator-(const typename AutoDiffBlock<Scalar>::V& lhs,
+                                    const AutoDiffBlock<Scalar>& rhs)
+    {
+        return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) - rhs;
+    }
+
+
+    /// Elementwise subtraction with constant on the right.
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator-(const AutoDiffBlock<Scalar>& lhs,
+                                    const typename AutoDiffBlock<Scalar>::V& rhs)
+    {
+        return lhs - AutoDiffBlock<Scalar>::constant(rhs, lhs.blockPattern());
+    }
+
+
+    /// Elementwise division with constant on the left.
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator/(const typename AutoDiffBlock<Scalar>::V& lhs,
+                                    const AutoDiffBlock<Scalar>& rhs)
+    {
+        return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) / rhs;
+    }
+
+
+    /// Elementwise division with constant on the right.
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator/(const AutoDiffBlock<Scalar>& lhs,
+                                    const typename AutoDiffBlock<Scalar>::V& rhs)
+    {
+        return lhs / AutoDiffBlock<Scalar>::constant(rhs, lhs.blockPattern());
+    }
+
+
+    /**
+     * @brief Operator for multiplication with a scalar on the right-hand side
+     *
+     * @param lhs The left-hand side AD forward block
+     * @param rhs The scalar to multiply with
+     * @return The product
+     */
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator*(const AutoDiffBlock<Scalar>& lhs,
+                                    const Scalar& rhs)
+    {
+        std::vector< typename AutoDiffBlock<Scalar>::M > jac;
+        jac.reserve( lhs.numBlocks() );
+        for (int block=0; block<lhs.numBlocks(); block++) {
+            jac.emplace_back( lhs.derivative()[block] * rhs );
+        }
+        return AutoDiffBlock<Scalar>::function( lhs.value() * rhs, jac );
+    }
+
+
+    /**
+     * @brief Operator for multiplication with a scalar on the left-hand side
+     *
+     * @param lhs The scalar to multiply with
+     * @param rhs The right-hand side AD forward block
+     * @return The product
+     */
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> operator*(const Scalar& lhs,
+                                    const AutoDiffBlock<Scalar>& rhs)
+    {
+        return rhs * lhs; // Commutative operation.
+    }
+
+
+} // namespace Opm
+
+
+
+#endif // OPM_AUTODIFFBLOCK_HEADER_INCLUDED

opm/polymer/fullyimplicit/AutoDiffHelpers.hpp

@@ -0,0 +1,513 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_AUTODIFFHELPERS_HEADER_INCLUDED
+#define OPM_AUTODIFFHELPERS_HEADER_INCLUDED
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+#include <iostream>
+#include <vector>
+
+namespace Opm
+{
+
+// -------------------- class HelperOps --------------------
+
+/// Contains vectors and sparse matrices that represent subsets or
+/// operations on (AD or regular) vectors of data.
+struct HelperOps
+{
+    typedef AutoDiffBlock<double>::M M;
+    typedef AutoDiffBlock<double>::V V;
+
+    /// A list of internal faces.
+    typedef Eigen::Array<int, Eigen::Dynamic, 1> IFaces;
+    IFaces internal_faces;
+
+    /// Extract for each internal face the difference of its adjacent cells' values (first - second).
+    M ngrad;
+    /// Extract for each face the difference of its adjacent cells' values (second - first).
+    M grad;
+    /// Extract for each face the average of its adjacent cells' values.
+    M caver;
+    /// Extract for each cell the sum of its adjacent interior faces' (signed) values.
+    M div;
+    /// Extract for each face the difference of its adjacent cells' values (first - second).
+    /// For boundary faces, one of the entries per row (corresponding to the outside) is zero.
+    M fullngrad;
+    /// Extract for each cell the sum of all its adjacent faces' (signed) values.
+    M fulldiv;
+
+    /// Constructs all helper vectors and matrices.
+    HelperOps(const UnstructuredGrid& grid)
+    {
+        const int nc = grid.number_of_cells;
+        const int nf = grid.number_of_faces;
+        // Define some neighbourhood-derived helper arrays.
+        typedef Eigen::Array<bool, Eigen::Dynamic, 1> OneColBool;
+        typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
+        typedef Eigen::Array<bool, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColBool;
+        TwoColInt nb = Eigen::Map<TwoColInt>(grid.face_cells, nf, 2);
+        // std::cout << "nb = \n" << nb << std::endl;
+        TwoColBool nbib = nb >= 0;
+        OneColBool ifaces = nbib.rowwise().all();
+        const int num_internal = ifaces.cast<int>().sum();
+        // std::cout << num_internal << " internal faces." << std::endl;
+        TwoColInt nbi(num_internal, 2);
+        internal_faces.resize(num_internal);
+        int fi = 0;
+        for (int f = 0; f < nf; ++f) {
+            if (ifaces[f]) {
+                internal_faces[fi] = f;
+                nbi.row(fi) = nb.row(f);
+                ++fi;
+            }
+        }
+        // std::cout << "nbi = \n" << nbi << std::endl;
+        // Create matrices.
+        ngrad.resize(num_internal, nc);
+        caver.resize(num_internal, nc);
+        typedef Eigen::Triplet<double> Tri;
+        std::vector<Tri> ngrad_tri;
+        std::vector<Tri> caver_tri;
+        ngrad_tri.reserve(2*num_internal);
+        caver_tri.reserve(2*num_internal);
+        for (int i = 0; i < num_internal; ++i) {
+            ngrad_tri.emplace_back(i, nbi(i,0), 1.0);
+            ngrad_tri.emplace_back(i, nbi(i,1), -1.0);
+            caver_tri.emplace_back(i, nbi(i,0), 0.5);
+            caver_tri.emplace_back(i, nbi(i,1), 0.5);
+        }
+        ngrad.setFromTriplets(ngrad_tri.begin(), ngrad_tri.end());
+        caver.setFromTriplets(caver_tri.begin(), caver_tri.end());
+        grad = -ngrad;
+        div = ngrad.transpose();
+        std::vector<Tri> fullngrad_tri;
+        fullngrad_tri.reserve(2*nf);
+        for (int i = 0; i < nf; ++i) {
+            if (nb(i,0) >= 0) {
+                fullngrad_tri.emplace_back(i, nb(i,0), 1.0);
+            }
+            if (nb(i,1) >= 0) {
+                fullngrad_tri.emplace_back(i, nb(i,1), -1.0);
+            }
+        }
+        fullngrad.resize(nf, nc);
+        fullngrad.setFromTriplets(fullngrad_tri.begin(), fullngrad_tri.end());
+        fulldiv = fullngrad.transpose();
+    }
+};
+
+// -------------------- upwinding helper class --------------------
+
+
+    /// Upwind selection in absence of counter-current flow (i.e.,
+    /// without effects of gravity and/or capillary pressure).
+    template <typename Scalar>
+    class UpwindSelector {
+    public:
+        typedef AutoDiffBlock<Scalar> ADB;
+
+        UpwindSelector(const UnstructuredGrid& g,
+                       const HelperOps&        h,
+                       const typename ADB::V&  ifaceflux)
+        {
+            typedef HelperOps::IFaces::Index IFIndex;
+            const IFIndex nif = h.internal_faces.size();
+            assert(nif == ifaceflux.size());
+
+            // Define selector structure.
+            typedef typename Eigen::Triplet<Scalar> Triplet;
+            std::vector<Triplet> s;  s.reserve(nif);
+            for (IFIndex iface = 0; iface < nif; ++iface) {
+                const int f  = h.internal_faces[iface];
+                const int c1 = g.face_cells[2*f + 0];
+                const int c2 = g.face_cells[2*f + 1];
+
+                assert ((c1 >= 0) && (c2 >= 0));
+
+                // Select upwind cell.
+                const int c = (ifaceflux[iface] >= 0) ? c1 : c2;
+
+                s.push_back(Triplet(iface, c, Scalar(1)));
+            }
+
+            // Assemble explicit selector operator.
+            select_.resize(nif, g.number_of_cells);
+            select_.setFromTriplets(s.begin(), s.end());
+        }
+
+        /// Apply selector to multiple per-cell quantities.
+        std::vector<ADB>
+        select(const std::vector<ADB>& xc) const
+        {
+            // Absence of counter-current flow means that the same
+            // selector applies to all quantities, 'x', defined per
+            // cell.
+            std::vector<ADB> xf;  xf.reserve(xc.size());
+            for (typename std::vector<ADB>::const_iterator
+                     b = xc.begin(), e = xc.end(); b != e; ++b)
+            {
+                xf.push_back(select_ * (*b));
+            }
+
+            return xf;
+        }
+
+        /// Apply selector to single per-cell ADB quantity.
+        ADB select(const ADB& xc) const
+        {
+            return select_*xc;
+        }
+
+        /// Apply selector to single per-cell constant quantity.
+        typename ADB::V select(const typename ADB::V& xc) const
+        {
+            return (select_*xc.matrix()).array();
+        }
+
+    private:
+        typename ADB::M select_;
+    };
+
+
+
+namespace {
+
+    template <typename Scalar, class IntVec>
+    Eigen::SparseMatrix<Scalar>
+    constructSubsetSparseMatrix(const int full_size, const IntVec& indices)
+    {
+        typedef Eigen::Triplet<Scalar> Tri;
+        const int subset_size = indices.size();
+        std::vector<Tri> triplets(subset_size);
+        for (int i = 0; i < subset_size; ++i) {
+            triplets[i] = Tri(i, indices[i], 1);
+        }
+        Eigen::SparseMatrix<Scalar> sub(subset_size, full_size);
+        sub.setFromTriplets(triplets.begin(), triplets.end());
+        return sub;
+    }
+
+    template <typename Scalar, class IntVec>
+    Eigen::SparseMatrix<Scalar>
+    constructSupersetSparseMatrix(const int full_size, const IntVec& indices)
+    {
+        return constructSubsetSparseMatrix<Scalar>(full_size, indices).transpose();
+    }
+
+} // anon namespace
+
+
+/// Returns x(indices).
+template <typename Scalar, class IntVec>
+AutoDiffBlock<Scalar>
+subset(const AutoDiffBlock<Scalar>& x,
+       const IntVec& indices)
+{
+    return constructSubsetSparseMatrix<Scalar>(x.value().size(), indices) * x;
+}
+
+
+
+/// Returns x(indices).
+template <typename Scalar, class IntVec>
+Eigen::Array<Scalar, Eigen::Dynamic, 1>
+subset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
+       const IntVec& indices)
+{
+    return (constructSubsetSparseMatrix<Scalar>(x.size(), indices) * x.matrix()).array();
+}
+
+
+
+/// Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
+template <typename Scalar, class IntVec>
+AutoDiffBlock<Scalar>
+superset(const AutoDiffBlock<Scalar>& x,
+         const IntVec& indices,
+         const int n)
+{
+    return constructSupersetSparseMatrix<Scalar>(n, indices) * x;
+}
+
+
+
+/// Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
+template <typename Scalar, class IntVec>
+Eigen::Array<Scalar, Eigen::Dynamic, 1>
+superset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
+         const IntVec& indices,
+         const int n)
+{
+    return constructSupersetSparseMatrix<Scalar>(n, indices) * x.matrix();
+}
+
+
+
+/// Construct square sparse matrix with the
+/// elements of d on the diagonal.
+/// Need to mark this as inline since it is defined in a header and not a template.
+inline
+AutoDiffBlock<double>::M
+spdiag(const AutoDiffBlock<double>::V& d)
+{
+    typedef AutoDiffBlock<double>::M M;
+
+    const int n = d.size();
+    M mat(n, n);
+    mat.reserve(Eigen::ArrayXi::Ones(n, 1));
+    for (M::Index i = 0; i < n; ++i) {
+        mat.insert(i, i) = d[i];
+    }
+
+    return mat;
+}
+
+
+
+
+    /// Selection. Choose first of two elements if selection basis element is nonnegative.
+    template <typename Scalar>
+    class Selector {
+    public:
+        typedef AutoDiffBlock<Scalar> ADB;
+
+        Selector(const typename ADB::V& selection_basis)
+        {
+            // Define selector structure.
+            const int n = selection_basis.size();
+            // Over-reserving so we do not have to count.
+            left_elems_.reserve(n);
+            right_elems_.reserve(n);
+            for (int i = 0; i < n; ++i) {
+                if (selection_basis[i] < 0.0) {
+                    right_elems_.push_back(i);
+                } else {
+                    left_elems_.push_back(i);
+                }
+            }
+        }
+
+        /// Apply selector to ADB quantities.
+        ADB select(const ADB& x1, const ADB& x2) const
+        {
+            if (right_elems_.empty()) {
+                return x1;
+            } else if (left_elems_.empty()) {
+                return x2;
+            } else {
+                return superset(subset(x1, left_elems_), left_elems_, x1.size())
+                    + superset(subset(x2, right_elems_), right_elems_, x2.size());
+            }
+        }
+
+        /// Apply selector to ADB quantities.
+        typename ADB::V select(const typename ADB::V& x1, const typename ADB::V& x2) const
+        {
+            if (right_elems_.empty()) {
+                return x1;
+            } else if (left_elems_.empty()) {
+                return x2;
+            } else {
+                return superset(subset(x1, left_elems_), left_elems_, x1.size())
+                    + superset(subset(x2, right_elems_), right_elems_, x2.size());
+            }
+        }
+
+    private:
+        std::vector<int> left_elems_;
+        std::vector<int> right_elems_;
+    };
+
+
+
+
+/// Returns the input expression, but with all Jacobians collapsed to one.
+inline
+AutoDiffBlock<double>
+collapseJacs(const AutoDiffBlock<double>& x)
+{
+    typedef AutoDiffBlock<double> ADB;
+    const int nb = x.numBlocks();
+    typedef Eigen::Triplet<double> Tri;
+    int nnz = 0;
+    for (int block = 0; block < nb; ++block) {
+        nnz += x.derivative()[block].nonZeros();
+    }
+    std::vector<Tri> t;
+    t.reserve(nnz);
+    int block_col_start = 0;
+    for (int block = 0; block < nb; ++block) {
+        const ADB::M& jac = x.derivative()[block];
+        for (ADB::M::Index k = 0; k < jac.outerSize(); ++k) {
+            for (ADB::M::InnerIterator i(jac, k); i ; ++i) {
+                t.push_back(Tri(i.row(),
+                                i.col() + block_col_start,
+                                i.value()));
+            }
+        }
+        block_col_start += jac.cols();
+    }
+    // Build final jacobian.
+    std::vector<ADB::M> jacs(1);
+    jacs[0].resize(x.size(), block_col_start);
+    jacs[0].setFromTriplets(t.begin(), t.end());
+    return ADB::function(x.value(), jacs);
+}
+
+
+
+
+/// Returns the vertical concatenation [ x; y ] of the inputs.
+inline
+AutoDiffBlock<double>
+vertcat(const AutoDiffBlock<double>& x,
+        const AutoDiffBlock<double>& y)
+{
+    const int nx = x.size();
+    const int ny = y.size();
+    const int n = nx + ny;
+    std::vector<int> xind(nx);
+    for (int i = 0; i < nx; ++i) {
+        xind[i] = i;
+    }
+    std::vector<int> yind(ny);
+    for (int i = 0; i < ny; ++i) {
+        yind[i] = nx + i;
+    }
+    return superset(x, xind, n) + superset(y, yind, n);
+}
+
+
+
+
+
+class Span
+{
+public:
+    explicit Span(const int num)
+    : num_(num),
+      stride_(1),
+      start_(0)
+    {
+    }
+    Span(const int num, const int stride, const int start)
+        : num_(num),
+          stride_(stride),
+          start_(start)
+    {
+    }
+    int operator[](const int i) const
+    {
+        assert(i >= 0 && i < num_);
+        return start_ + i*stride_;
+    }
+    int size() const
+    {
+        return num_;
+    }
+
+
+    class SpanIterator
+    {
+    public:
+        SpanIterator(const Span* span, const int index)
+            : span_(span),
+              index_(index)
+        {
+        }
+        SpanIterator operator++()
+        {
+            ++index_;
+            return *this;
+        }
+        SpanIterator operator++(int)
+        {
+            SpanIterator before_increment(*this);
+            ++index_;
+            return before_increment;
+        }
+        bool operator<(const SpanIterator& rhs) const
+        {
+            assert(span_ == rhs.span_);
+            return index_ < rhs.index_;
+        }
+        bool operator==(const SpanIterator& rhs) const
+        {
+            assert(span_ == rhs.span_);
+            return index_ == rhs.index_;
+        }
+        bool operator!=(const SpanIterator& rhs) const
+        {
+            assert(span_ == rhs.span_);
+            return index_ != rhs.index_;
+        }
+        int operator*()
+        {
+            return (*span_)[index_];
+        }
+    private:
+        const Span* span_;
+        int index_;
+    };
+
+    typedef SpanIterator iterator;
+    typedef SpanIterator const_iterator;
+
+    SpanIterator begin() const
+    {
+        return SpanIterator(this, 0);
+    }
+
+    SpanIterator end() const
+    {
+        return SpanIterator(this, num_);
+    }
+
+    bool operator==(const Span& rhs)
+    {
+        return num_ == rhs.num_ && start_ == rhs.start_ && stride_ == rhs.stride_;
+    }
+
+private:
+    const int num_;
+    const int stride_;
+    const int start_;
+};
+
+
+
+/// Return a vector of (-1.0, 0.0 or 1.0), depending on sign per element.
+inline Eigen::ArrayXd sign (const Eigen::ArrayXd& x)
+{
+    const int n = x.size();
+    Eigen::ArrayXd retval(n);
+    for (int i = 0; i < n; ++i) {
+        retval[i] = x[i] < 0.0 ? -1.0 : (x[i] > 0.0 ? 1.0 : 0.0);
+    }
+    return retval;
+}
+
+} // namespace Opm
+
+#endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED

opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp

@@ -0,0 +1,449 @@
+/**/
+
+#include <opm/autodiff/polymer/FullyImplicitTwoPhaseSolver.hpp>
+
+#include <opm/core/pressure/tpfa/trans_tpfa.h>
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/polymer/IncompPropsAdInterface.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include <iomanip>
+#include <Eigen/Eigen>
+#include <algorithm>
+namespace Opm {
+
+namespace {
+    
+    std::vector<int>
+    buildAllCells(const int nc)
+    {
+        std::vector<int> all_cells(nc);
+        for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
+
+        return all_cells;
+    }
+    struct Chop01 {
+        double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); }
+    };
+
+}//anonymous namespace
+
+
+
+
+
+typedef AutoDiffBlock<double> ADB;
+typedef ADB::V V;
+typedef ADB::M M;
+typedef Eigen::Array<double,
+                     Eigen::Dynamic,
+                     Eigen::Dynamic,
+                     Eigen::RowMajor> DataBlock;
+
+
+
+
+
+    FullyImplicitTwoPhaseSolver::
+    FullyImplicitTwoPhaseSolver(const UnstructuredGrid&         grid,
+                                const IncompPropsAdInterface&   fluid,
+                                const LinearSolverInterface&    linsolver)
+        : grid_ (grid)
+        , fluid_(fluid)
+        , linsolver_(linsolver)
+        , cells_ (buildAllCells(grid.number_of_cells))
+        , ops_(grid)
+        , residual_(std::vector<ADB>(fluid.numPhases(), ADB::null()))
+     {
+     }
+
+
+
+
+
+    void
+    FullyImplicitTwoPhaseSolver::
+    step(const double   dt,
+         TwophaseState& x,
+         const std::vector<double>& src)
+    {
+        
+        V pvol(grid_.number_of_cells);
+        // Pore volume
+        const typename V::Index nc = grid_.number_of_cells;
+        V rho = V::Constant(pvol.size(), 1, *fluid_.porosity());
+        std::transform(grid_.cell_volumes, grid_.cell_volumes + nc,
+                       rho.data(), pvol.data(),
+                       std::multiplies<double>());
+
+        const V pvdt = pvol / dt;
+
+        const SolutionState old_state = constantState(x);
+        const double atol  = 1.0e-12;
+        const double rtol  = 5.0e-8;
+        const int    maxit = 15;
+
+        assemble(pvdt, old_state, x, src);
+
+        const double r0  = residualNorm();
+        int          it  = 0;
+        std::cout << "\nIteration         Residual\n"
+                  << std::setw(9) << it << std::setprecision(9)
+                  << std::setw(18) << r0 << std::endl;
+        bool resTooLarge = r0 > atol;
+        while (resTooLarge && (it < maxit)) {
+            const V dx = solveJacobianSystem();
+            updateState(dx, x);
+
+            assemble(pvdt, old_state, x, src);
+
+            const double r = residualNorm();
+
+            resTooLarge = (r > atol) && (r > rtol*r0);
+
+            it += 1;
+            std::cout << std::setw(9) << it << std::setprecision(9)
+                      << std::setw(18) << r << std::endl;
+        }
+
+        if (resTooLarge) {
+            std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
+            // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
+        }
+    }
+
+
+
+
+
+    FullyImplicitTwoPhaseSolver::SolutionState::SolutionState(const int np)
+        : pressure   (    ADB::null())
+        , saturation (np, ADB::null())
+    {
+    }
+
+
+
+
+
+    FullyImplicitTwoPhaseSolver::SolutionState
+    FullyImplicitTwoPhaseSolver::constantState(const TwophaseState& x)
+    {
+        const int nc = grid_.number_of_cells;
+        const int np = x.numPhases();
+
+        SolutionState state(np);
+
+        // Pressure.
+        assert (not x.pressure().empty());
+        const V p = Eigen::Map<const V>(& x.pressure()[0], nc);
+        state.pressure = ADB::constant(p);
+
+        // Saturation.
+        assert (not x.saturation().empty());
+        const DataBlock s_all = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
+        for (int phase = 0; phase < np; ++phase) {
+            state.saturation[phase] = ADB::constant(s_all.col(phase));
+   //     state.saturation[1] = ADB::constant(s_all.col(1));
+        }
+        return state;
+    }
+
+
+
+
+
+    FullyImplicitTwoPhaseSolver::SolutionState
+    FullyImplicitTwoPhaseSolver::variableState(const TwophaseState& x)
+    {
+        const int nc = grid_.number_of_cells;
+        const int np = x.numPhases();
+
+        std::vector<V> vars0;
+        vars0.reserve(np); 
+
+        // Initial pressure.
+        assert (not x.pressure().empty());
+        const V p = Eigen::Map<const V>(& x.pressure()[0], nc);
+        vars0.push_back(p);
+
+        // Initial saturation.
+        assert (not x.saturation().empty());
+        const DataBlock s_all = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
+        const V sw = s_all.col(0);
+        vars0.push_back(sw);
+
+
+        std::vector<ADB> vars = ADB::variables(vars0);
+
+        SolutionState state(np);
+
+        // Pressure.
+        int nextvar = 0;
+        state.pressure = vars[ nextvar++ ];
+
+        // Saturation.
+        const std::vector<int>& bpat = vars[0].blockPattern();
+        {
+            ADB so  = ADB::constant(V::Ones(nc, 1), bpat);
+            ADB& sw = vars[ nextvar++ ];
+            state.saturation[0] = sw;
+            so = so - sw;
+            state.saturation[1] = so;
+        }
+
+        assert(nextvar == int(vars.size()));
+
+        return state;
+    }
+  
+
+
+
+
+    void
+    FullyImplicitTwoPhaseSolver::
+    assemble(const V&             pvdt,
+             const SolutionState& old_state,
+             const TwophaseState& x   ,
+             const std::vector<double>&  src)
+    {
+        // Create the primary variables.
+        const SolutionState state = variableState(x);
+
+        // -------- Mass balance equations --------
+        const V trans = subset(transmissibility(), ops_.internal_faces);
+        const std::vector<ADB> kr = computeRelPerm(state);
+        for (int phase = 0; phase < fluid_.numPhases(); ++phase) {
+            const ADB mflux = computeMassFlux(phase, trans, kr, state);
+            ADB source = accumSource(phase, kr, src);
+            residual_[phase] =
+                pvdt*(state.saturation[phase] - old_state.saturation[phase])
+                + ops_.div*mflux - source;
+        }
+
+    }
+   
+
+
+
+
+
+    ADB 
+    FullyImplicitTwoPhaseSolver::accumSource(const int phase,
+                                                 const std::vector<ADB>& kr,
+                                                 const std::vector<double>& src) const
+    {
+        //extract the source to out and in source.
+        std::vector<double> outsrc;
+        std::vector<double> insrc;
+        std::vector<double>::const_iterator it;
+        for (it = src.begin(); it != src.end(); ++it) {
+            if (*it < 0) {
+                outsrc.push_back(*it);
+                insrc.push_back(0.0);
+            } else if (*it > 0) {
+                insrc.push_back(*it);
+                outsrc.push_back(0.0);
+            } else {
+                outsrc.emplace_back(0);
+                insrc.emplace_back(0);
+            }
+        }
+        const V source = Eigen::Map<const V>(& src[0], grid_.number_of_cells);
+        const V outSrc = Eigen::Map<const V>(& outsrc[0], grid_.number_of_cells);
+        const V inSrc = Eigen::Map<const V>(& insrc[0], grid_.number_of_cells);
+        
+        // compute the out-fracflow.
+        ADB f_out = computeFracFlow(phase, kr);
+        // compute the in-fracflow.
+        V f_in;
+        if (phase == 1) {
+            f_in = V::Zero(grid_.number_of_cells);
+        } else if (phase == 0) {
+            f_in = V::Ones(grid_.number_of_cells);
+        }
+        return f_out * outSrc + f_in * inSrc;
+     }
+
+
+
+
+
+    ADB
+    FullyImplicitTwoPhaseSolver::computeFracFlow(int                     phase,
+                                                 const std::vector<ADB>& kr) const
+    {
+        const double* mus = fluid_.viscosity();
+        ADB  mob_phase = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]);
+        ADB  mob_wat = kr[0] / V::Constant(kr[0].size(), 1, mus[0]);
+        ADB  mob_oil= kr[1] / V::Constant(kr[1].size(), 1, mus[1]);
+        ADB  total_mob = mob_wat + mob_oil;
+        ADB f = mob_phase / total_mob;
+
+        return f;
+    }
+
+
+
+
+
+    V 
+    FullyImplicitTwoPhaseSolver::solveJacobianSystem() const
+    {
+        const int np = fluid_.numPhases();
+    	if (np != 2) {
+	        OPM_THROW(std::logic_error, "Only two-phase ok in FullyImplicitTwoPhaseSolver.");
+	    }
+	    ADB mass_res = collapseJacs(vertcat(residual_[0], residual_[1]));
+        const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = mass_res.derivative()[0];
+        V dx(V::Zero(mass_res.size()));
+        Opm::LinearSolverInterface::LinearSolverReport rep
+            = linsolver_.solve(matr.rows(), matr.nonZeros(),
+                               matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
+                               mass_res.value().data(), dx.data());
+        if (!rep.converged) {
+            OPM_THROW(std::runtime_error,
+                      "FullyImplicitBlackoilSolver::solveJacobianSystem(): "
+                      "Linear solver convergence failure.");
+        }
+        return dx;
+    }
+
+
+
+
+
+    void FullyImplicitTwoPhaseSolver::updateState(const V& dx,
+                                                  TwophaseState& state) const
+    {
+        const int np = fluid_.numPhases();
+        const int nc = grid_.number_of_cells;
+        const V null;
+        assert(null.size() == 0);
+        const V zero = V::Zero(nc);
+        const V one = V::Constant(nc, 1.0);
+
+        // Extract parts of dx corresponding to each part.
+        const V dp = subset(dx, Span(nc));
+        int varstart = nc;
+        const V dsw = subset(dx, Span(nc, 1, varstart));
+        varstart += dsw.size();
+        assert(varstart == dx.size());
+
+        // Pressure update.
+        const V p_old = Eigen::Map<const V>(&state.pressure()[0], nc);
+        const V p = p_old - dp;
+        std::copy(&p[0], &p[0] + nc, state.pressure().begin());
+
+
+        // Saturation updates.
+        const double dsmax = 0.3;
+        const DataBlock s_old = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
+        V so = one;
+        const V sw_old = s_old.col(0);
+        const V dsw_limited = sign(dsw) * dsw.abs().min(dsmax);
+        const V sw = (sw_old - dsw_limited).unaryExpr(Chop01());
+        so -= sw;
+        for (int c = 0; c < nc; ++c) {
+            state.saturation()[c*np] = sw[c];
+        }
+        for (int c = 0; c < nc; ++c) {
+            state.saturation()[c*np + 1] = so[c];
+        }
+
+    }
+   
+
+
+
+
+    std::vector<ADB>
+    FullyImplicitTwoPhaseSolver::computeRelPerm(const SolutionState& state) const
+    {
+
+        const ADB sw = state.saturation[0];
+        const ADB so = state.saturation[1];
+
+        return fluid_.relperm(sw, so, cells_);
+    }
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    ADB
+    FullyImplicitTwoPhaseSolver::computeMassFlux(const int               phase ,
+                                                 const V&                trans,
+                                                 const std::vector<ADB>& kr    ,
+                                                 const SolutionState&    state ) const
+    {
+//        const ADB tr_mult = transMult(state.pressure);
+        const double* mus = fluid_.viscosity();
+        ADB  mob = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]);
+
+        const ADB dp = ops_.ngrad * state.pressure;
+        const ADB head = trans * dp;
+
+        UpwindSelector<double> upwind(grid_, ops_, head.value());
+
+        return upwind.select(mob) * head;
+    }
+  
+
+   
+   
+   
+    double
+    FullyImplicitTwoPhaseSolver::residualNorm() const
+    {
+        double r = 0;
+        for (std::vector<ADB>::const_iterator
+                 b = residual_.begin(),
+                 e = residual_.end();
+             b != e; ++b)
+        {
+            r = std::max(r, (*b).value().matrix().norm());
+        }
+
+        return r;
+    }
+   
+   
+   
+   
+   
+    V
+    FullyImplicitTwoPhaseSolver::transmissibility() const
+    {
+        const V::Index nc = grid_.number_of_cells;
+        V htrans(grid_.cell_facepos[nc]);
+        V trans(grid_.cell_facepos[nc]);
+        UnstructuredGrid* ug = const_cast<UnstructuredGrid*>(& grid_);
+        tpfa_htrans_compute(ug, fluid_.permeability(), htrans.data());
+        tpfa_trans_compute (ug, htrans.data()     , trans.data());
+        
+        return trans;
+    }
+
+
+
+
+
+
+}//namespace Opm

opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.hpp

@@ -0,0 +1,91 @@
+/**/
+
+#ifndef OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED
+#define OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/polymer/IncompPropsAdInterface.hpp>
+#include <opm/core/pressure/tpfa/trans_tpfa.h>
+
+
+struct UnstructuredGrid;
+namespace Opm {
+//    struct HelperOps;
+    class LinearSolverInterface;
+    class TwophaseState;
+
+    
+    class FullyImplicitTwoPhaseSolver
+    {
+    public:
+        FullyImplicitTwoPhaseSolver(const UnstructuredGrid&        grid,
+                                    const IncompPropsAdInterface&  fluid,
+                                    const LinearSolverInterface&    linsolver);
+
+        void step(const double   dt,
+                  TwophaseState& state,
+                  const std::vector<double>& src);
+    private:
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        typedef ADB::M M;
+        typedef Eigen::Array<double,
+                             Eigen::Dynamic,
+                             Eigen::Dynamic,
+                             Eigen::RowMajor> DataBlock;
+        struct SolutionState {
+            SolutionState(const int np);
+            ADB             pressure;
+            std::vector<ADB> saturation;
+        };
+        const UnstructuredGrid&         grid_;
+        const IncompPropsAdInterface&   fluid_;
+        const LinearSolverInterface&    linsolver_;
+        const std::vector<int>          cells_;
+        HelperOps                       ops_;
+        std::vector<ADB>                residual_;
+       
+
+        SolutionState
+        constantState(const TwophaseState& x);
+        SolutionState
+        variableState(const TwophaseState& x);
+        void
+        assemble(const V&               pvdt,
+                 const SolutionState&   old_state,
+                 const TwophaseState&  x,
+                 const std::vector<double>& src);
+        V solveJacobianSystem() const;
+        void updateState(const V&             dx,
+                         TwophaseState& x) const;
+        std::vector<ADB>
+        computeRelPerm(const SolutionState& state) const;
+        V
+        transmissibility() const;
+        ADB
+        computeFracFlow(int    phase,
+                        const std::vector<ADB>& kr) const;
+        ADB 
+        accumSource(const int phase,
+                    const std::vector<ADB>& kr,
+                    const std::vector<double>& src) const;
+        ADB
+        computeMassFlux(const int               phase,
+                        const V&                trans,
+                        const std::vector<ADB>& kr,
+                        const SolutionState&    state) const;
+        double
+        residualNorm() const;
+
+        ADB
+        rockPorosity(const ADB& p) const;
+        ADB
+        rockPermeability(const ADB& p) const;
+        const double
+        fluidDensity(const int phase) const;
+        ADB
+        transMult(const ADB& p) const;
+    };
+} // namespace Opm
+#endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED

opm/polymer/fullyimplicit/IncompPropsAd.cpp

@@ -0,0 +1,5 @@
+/*
+Author : Liu Ming 
+Date   : 2013-11-28 in Oslo
+Email  : miliu@statoil.com
+*/

opm/polymer/fullyimplicit/IncompPropsAd.hpp

@@ -0,0 +1,59 @@
+/*
+Author : Liu Ming
+Data   : 2013-11-28  in Oslo
+Email  : miliu@statoil.com
+
+Properties for incompressible immiscible two-phase flow
+*/
+
+#ifndef OPM_INCOMPROPSAD_HEADER_INCLUDED
+#define  OPM_INCOMPROPSAD_HEADER_INCLUDED
+#include <opm/core/props/IncompPropertiesBasic.hpp>
+#include <opm/core/props/rock/RockBasic.hpp>
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/core/props/pvt/PvtPropertiesBasic.hpp>
+#include <opm/core/props/satfunc/SaturationPropsBasic.hpp>
+
+
+namespace Opm
+{
+
+//    class BlackoilPhases;
+
+    class IncompropsAd : public IncompPropertiesBasic
+    {
+    public:
+        IncomPropsAd(const parameter::ParameterGroup& param,
+                     const int dim,
+                     const int num_cells);
+        IncompPropsAd(const int num_phases,
+                      const SaturationPropsBasic::RelPermFunc& relpermfunc,
+                      const std::vector<double>& rho,
+                      const std::vector<double>& mu,
+                      const double porosity,
+                      const double permeability,
+                      const int dim,
+                      const int num_cells);
+
+        ~IncompPropsAd();
+
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        V relperm(const V& sw,
+                  const V& so,
+                  const std::vector<int>& cells);
+        ADB relperm(const ADB& sw,
+                    const ADB& so,
+                    const std::vector<int>& cells);
+
+        V capPress(const V& sw,
+                   const V& so,
+                   const std::vector<int>& cells);
+        ADB capPress(const ADB& sw,
+                     const ADB& so,
+                     const std::vector<int>& cells);
+
+    }
+}
+
+#endif//  OPM_INCOMPROPSAD_HEADER_INCLUDED

opm/polymer/fullyimplicit/IncompPropsAdBasic.cpp

@@ -0,0 +1,172 @@
+/**/
+
+#include <opm/autodiff/polymer/IncompPropsAdBasic.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <iostream>
+
+namespace Opm
+{
+    IncompPropsAdBasic::IncompPropsAdBasic(const parameter::ParameterGroup& param,
+                                           const int dim,
+                                           const int num_cells)
+    {
+        double poro = param.getDefault("porosity", 1.0);
+        using namespace Opm::unit;
+        using namespace Opm::prefix;
+        double perm = param.getDefault("permeability", 100) * milli * darcy;
+        rock_.init(dim, num_cells, poro, perm);
+        pvt_.init(param);
+        satprops_.init(param);
+        if (pvt_.numPhases() != satprops_.numPhases()) {
+            OPM_THROW(std::runtime_error, "IncompPropsAdBasic::IncompPropsAdBasic() - Inconsistent number of phases in pvt data ("
+                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
+        }
+        viscosity_.resize(pvt_.numPhases());
+        pvt_.mu(1, 0, 0, &viscosity_[0]);
+    }
+
+
+    IncompPropsAdBasic::IncompPropsAdBasic(const int num_phases,
+                                                 const  SaturationPropsBasic::RelPermFunc& relpermfunc,
+                                                 const std::vector<double>&  rho,
+                                                 const std::vector<double>& mu,
+                                                 const double por, //porosity
+                                                 const double perm,
+                                                 const int dim,
+                                                 const int num_cells)
+    {
+        rock_.init(dim, num_cells, por, perm);
+        pvt_.init(num_phases, rho, mu);
+        satprops_.init(num_phases, relpermfunc);
+        if (pvt_.numPhases() != satprops_.numPhases()) {
+            OPM_THROW(std::runtime_error, "IncompPropsAdBasic::IncompPropsAdBasic() - Inconsistent number of phases in pvt data ("
+                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
+        }
+        viscosity_.resize(pvt_.numPhases());
+        pvt_.mu(1, 0, 0, &viscosity_[0]);
+    }
+    IncompPropsAdBasic::~IncompPropsAdBasic()
+    {
+    }
+
+    /// \return   D, the number of spatial dimensions.
+    int IncompPropsAdBasic::numDimensions() const
+    {
+        return rock_.numDimensions();
+    }
+
+    /// \return   N, the number of cells.
+    int IncompPropsAdBasic::numCells() const
+    {
+        return rock_.numCells();
+    }
+
+    /// \return   Array of N porosity values.
+    const double* IncompPropsAdBasic::porosity() const
+    {
+        return rock_.porosity();
+    }
+
+    /// \return   Array of ND^2 permeability values.
+    ///           The D^2 permeability values for a cell are organized as a matrix,
+    ///           which is symmetric (so ordering does not matter).
+    const double* IncompPropsAdBasic::permeability() const
+    {
+        return rock_.permeability();
+    }
+
+
+    // ---- Fluid interface ----
+
+    /// \return   P, the number of phases (also the number of components).
+    int IncompPropsAdBasic::numPhases() const
+    {
+        return pvt_.numPhases();
+    }
+
+    /// \return Array of P viscosity values.
+    const double* IncompPropsAdBasic::viscosity() const
+    {
+        return &viscosity_[0];
+    }
+    /// \return Array of P density values.
+    const double* IncompPropsAdBasic::density() const
+    {
+        return pvt_.surfaceDensities();
+    }
+
+    const double* IncompPropsAdBasic::surfaceDensity() const
+    {
+        return pvt_.surfaceDensities();
+    }
+    typedef IncompPropsAdBasic::ADB ADB;
+    typedef IncompPropsAdBasic::V V;
+    typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Block;
+    typedef std::vector<int> Cells;
+
+    std::vector<V> 
+    IncompPropsAdBasic::relperm(const V& sw,
+                                const V& so,
+                                const Cells& cells) const
+    {
+        const int n = cells.size();
+        const int np = numPhases();
+        Block s_all(n, np);
+        assert(sw.size() == n && so.size() == n);
+        s_all.col(0) = sw;
+        s_all.col(1) = so;
+        Block kr(n, np);
+//        satprops_.relperm(n, s_all.data(), cells.data(), kr.data(), 0);
+        satprops_.relperm(n, s_all.data(), kr.data(), 0);
+
+        std::vector<V> relperms;
+        relperms.reserve(2);
+        for (int phase = 0; phase < 2; ++phase) {
+                relperms.emplace_back(kr.col(phase));
+        }
+        return relperms;
+    }
+
+    std::vector<ADB>    
+    IncompPropsAdBasic::relperm(const ADB& sw,
+                                const ADB& so,
+                                const Cells& cells) const
+    {
+        const int n = cells.size();
+        const int np = numPhases();
+        Block s_all(n, np);
+        assert(sw.size() == n && so.size() == n);
+        s_all.col(0) = sw.value();
+        s_all.col(1) = so.value();
+        Block kr(n, np);
+        Block dkr(n, np*np);
+//        satprops_.relperm(n, s_all.data(), cells.data(), kr.data(), dkr.data());
+        satprops_.relperm(n, s_all.data(), kr.data(), dkr.data());
+        const int num_blocks = so.numBlocks();
+        std::vector<ADB> relperms;
+        relperms.reserve(2);
+        typedef const ADB* ADBPtr;
+        ADBPtr s[2] = { &sw, &so };
+        for (int phase1 = 0; phase1 < 2; ++phase1) {
+            const int phase1_pos = phase1;
+            std::vector<ADB::M> jacs(num_blocks);
+            for (int block = 0; block < num_blocks; ++block) {
+                jacs[block] = ADB::M(n, s[phase1]->derivative()[block].cols());
+            }
+            for (int phase2 = 0; phase2 < 2; ++phase2) {
+                const int phase2_pos = phase2;
+                // Assemble dkr1/ds2.
+                const int column = phase1_pos + np*phase2_pos; // Recall: Fortran ordering from props_.relperm()
+                ADB::M dkr1_ds2_diag = spdiag(dkr.col(column));
+                for (int block = 0; block < num_blocks; ++block) {
+                    jacs[block] += dkr1_ds2_diag * s[phase2]->derivative()[block];
+                }
+            }
+                relperms.emplace_back(ADB::function(kr.col(phase1_pos), jacs));
+        }
+        return relperms;
+    }
+}

opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp

@@ -0,0 +1,54 @@
+/**/
+
+#ifndef OPM_INCOMPPROPSADBASIC_HEADER_INCLUDED
+#define OPM_INCOMPPROPSADBASIC_HEADER_INCLUDED
+
+#include <opm/autodiff/polymer/IncompPropsAdInterface.hpp>
+#include <opm/core/props/rock/RockBasic.hpp>
+#include <opm/core/props/pvt/PvtPropertiesBasic.hpp>
+#include <opm/core/props/satfunc/SaturationPropsBasic.hpp>
+#include <opm/autodiff/AutoDiffBlock.hpp>
+namespace Opm
+{
+    class IncompPropsAdBasic : public IncompPropsAdInterface
+    {
+    public:
+        IncompPropsAdBasic(const parameter::ParameterGroup& param,
+                           const int dim,
+                           const int num_cells);
+        IncompPropsAdBasic(const int num_phases,
+                              const SaturationPropsBasic::RelPermFunc& relpermfunc,
+                              const std::vector<double>& rho,
+                              const std::vector<double>& mu,
+                              const double porosity,
+                              const double permeability,
+                              const int dim,
+                              const int num_cells);
+
+         ~IncompPropsAdBasic();
+         int numDimensions() const;
+         int numCells() const;
+         const double* porosity() const;
+         const double* permeability() const;
+         
+         int numPhases() const;
+         const double* viscosity() const;
+         const double* density() const;
+         const double* surfaceDensity() const;
+
+         typedef AutoDiffBlock<double> ADB;
+         typedef ADB::V V;
+         std::vector<V> relperm(const V& sw,
+                                const V& so,
+                                const std::vector<int>& cells) const;
+         std::vector<ADB> relperm(const ADB& sw,
+                     const ADB& so,
+                     const std::vector<int>& cells) const;
+    private:
+        RockBasic rock_;
+        PvtPropertiesBasic pvt_;
+        SaturationPropsBasic satprops_;
+        std::vector<double> viscosity_;
+    };
+}
+#endif // OPM_INCOMPPROPSADBASIC_HEADER_INCLUDED

opm/polymer/fullyimplicit/IncompPropsAdFromDeck.cpp

@@ -0,0 +1,136 @@
+/**/
+#include <opm/autodiff/polymer/IncompPropsAdFromDeck.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <iostream>
+
+namespace Opm
+{
+    IncompPropsAdFromDeck::IncompPropsAdFromDeck(const EclipseGridParser& deck,
+                                                 const UnstructuredGrid& grid)
+    {
+        rock_.init(deck, grid);
+        pvt_.init(deck);
+        satprops_.init(deck, grid, 200);
+        if (pvt_.numPhases() != satprops_.numPhases()) {
+            OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "
+                  "Inconsistent number of phases in pvt data (" << pvt_.numPhases()
+                  << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
+        }
+    }
+    IncompPropsAdFromDeck::~IncompPropsAdFromDeck()
+    {
+    }
+    // rock interface
+    int IncompPropsAdFromDeck::numDimensions() const
+    {
+        return rock_.numDimensions();
+    }
+    int IncompPropsAdFromDeck::numCells() const
+    {
+        return rock_.numCells();
+    }
+
+    const double* IncompPropsAdFromDeck::porosity() const
+    {
+        return rock_.porosity();
+    }
+    const double* IncompPropsAdFromDeck::permeability() const
+    {
+        return rock_.permeability();
+    }
+
+    // fluid interface
+    int IncompPropsAdFromDeck::numPhases() const
+    {
+        return pvt_.numPhases();
+    }
+
+
+    const double* IncompPropsAdFromDeck::viscosity() const
+    {
+        return pvt_.viscosity();
+    }
+
+
+    const double* IncompPropsAdFromDeck::density() const
+    {
+        return pvt_.reservoirDensities();
+    }
+
+
+    const double* IncompPropsAdFromDeck::surfaceDensity() const
+    {
+        return pvt_.surfaceDensities();
+    }
+
+
+    typedef IncompPropsAdFromDeck::ADB ADB;
+    typedef IncompPropsAdFromDeck::V V;
+    typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Block;
+
+
+    std::vector<V> 
+    IncompPropsAdFromDeck::relperm(const V& sw,
+                                   const V& so,
+                                   const Cells& cells) const
+    {
+        const int n = cells.size();
+        const int np = numPhases();
+        Block s_all(n, np);
+        assert(sw.size() == n && so.size() == n);
+        s_all.col(0) = sw;
+        s_all.col(1) = so;
+        Block kr(n, np);
+        satprops_.relperm(n, s_all.data(), cells.data(), kr.data(), 0);
+        std::vector<V> relperms;
+        relperms.reserve(np);
+        for (int phase = 0; phase < np; ++phase) {
+            relperms.emplace_back(kr.col(phase));
+        }
+        return relperms;
+    }
+
+
+
+    std::vector<ADB> 
+    IncompPropsAdFromDeck::relperm(const ADB& sw,
+                                   const ADB& so,
+                                   const Cells& cells) const
+    {
+        const int n = cells.size();
+        const int np = numPhases();
+        Block s_all(n, np);
+        assert(sw.size() == n && so.size() == n);
+        s_all.col(0) = sw.value();
+        s_all.col(1) = so.value();
+        Block kr(n, np);
+        Block dkr(n, np*np);
+        satprops_.relperm(n, s_all.data(), cells.data(), kr.data(), dkr.data());
+        const int num_blocks = so.numBlocks();
+        std::vector<ADB> relperms;
+        relperms.reserve(np);
+        typedef const ADB* ADBPtr;
+        ADBPtr s[2] = { &sw, &so };
+        for (int phase1 = 0; phase1 < np; ++phase1) {
+            const int phase1_pos = phase1;
+            std::vector<ADB::M> jacs(num_blocks);
+            for (int block = 0; block < num_blocks; ++block) {
+                jacs[block] = ADB::M(n, s[phase1]->derivative()[block].cols());
+            }
+            for (int phase2 = 0; phase2 < np; ++phase2) {
+                const int phase2_pos = phase2;
+                // Assemble dkr1/ds2.
+                const int column = phase1_pos + np*phase2_pos; // Recall: Fortran ordering from props_.relperm()
+                ADB::M dkr1_ds2_diag = spdiag(dkr.col(column));
+                for (int block = 0; block < num_blocks; ++block) {
+                    jacs[block] += dkr1_ds2_diag * s[phase2]->derivative()[block];
+                }
+            }
+                relperms.emplace_back(ADB::function(kr.col(phase1_pos), jacs));
+        }
+        return relperms;
+    }
+} //namespace Opm
+

opm/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp

@@ -0,0 +1,57 @@
+/**/
+#ifndef OPM_INCOMPPROPSADFROMDECK_HEADER_INCLUDED
+#define OPM_INCOMPPROPSADFROMDECK_HEADER_INCLUDED
+#include <opm/autodiff/polymer/IncompPropsAdInterface.hpp>
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/core/props/rock/RockFromDeck.hpp>
+#include <opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp>
+#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
+#include <opm/core/io/eclipse/EclipseGridParser.hpp>
+
+#include <boost/shared_ptr.hpp>
+#include <boost/scoped_ptr.hpp>
+
+struct UnstructuredGrid;
+namespace Opm
+{
+    class IncompPropsAdFromDeck : public IncompPropsAdInterface
+    {
+    public:
+        IncompPropsAdFromDeck(const EclipseGridParser& deck,
+                              const UnstructuredGrid&  grid);  
+        ~IncompPropsAdFromDeck();
+
+        //--Rock interface--
+        int numDimensions() const;
+        int numCells() const;
+        const double* porosity() const;
+        const double* permeability() const;
+
+        // --Fluid interface--
+        int numPhases() const;
+        const double* viscosity() const;
+        const double* density() const;
+        const double* surfaceDensity() const;
+
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        typedef std::vector<int> Cells;
+
+
+        std::vector<V> relperm(const V& sw,
+                               const V& so,
+                               const Cells& cells) const;
+
+        std::vector<ADB> relperm(const ADB& sw,
+                                 const ADB& so,
+                                 const Cells& cells) const;
+
+    private:
+        RockFromDeck rock_;
+        PvtPropertiesIncompFromDeck pvt_;
+        SaturationPropsFromDeck<SatFuncStone2Uniform> satprops_;
+    };
+
+} //namespace Opm
+
+#endif// OPM_INCOMPPROPSADFROMDECK_HEADER_INCLUDED

opm/polymer/fullyimplicit/IncompPropsAdInterface.cpp

@@ -0,0 +1,6 @@
+
+#include <opm/autodiff/polymer/IncompPropsAdInterface.hpp>
+
+Opm::IncompPropsAdInterface::~IncompPropsAdInterface()
+{
+}

opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp

@@ -0,0 +1,89 @@
+/*
+    
+*/
+
+#ifndef OPM_INCOMPPROPSADINTERFACE_HEADER_INCLUDED
+#define OPM_INCOMPPROPSADINTERFACE_HEADER_INCLUDED
+#include <opm/autodiff/AutoDiffBlock.hpp>
+
+namespace Opm
+{
+    class IncompPropsAdInterface
+    {
+    public:
+        virtual ~IncompPropsAdInterface();
+
+        ////////////////////////////
+        //      Rock interface    //
+        ////////////////////////////
+
+        /// \return   D, the number of spatial dimensions.
+        virtual int numDimensions() const = 0;
+
+        /// \return   N, the number of cells.
+        virtual int numCells() const = 0;
+
+        /// \return   Array of N porosity values.
+        virtual const double* porosity() const = 0;
+
+        /// \return   Array of ND^2 permeability values.
+        ///           The D^2 permeability values for a cell are organized as a matrix,
+        ///           which is symmetric (so ordering does not matter).
+        virtual const double* permeability() const = 0;
+
+
+        ////////////////////////////
+        //      Fluid interface   //
+        ////////////////////////////
+
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        typedef ADB::M M;
+        typedef std::vector<int> Cells;
+
+        /// \return   Number of active phases (also the number of components).
+        virtual int numPhases() const = 0;
+
+        // ------ Density ------
+
+        /// Densities of stock components at surface conditions.
+        /// \return Array of 2 density values.
+        virtual const double* surfaceDensity() const = 0;
+    
+
+        // ------ Viscosity ------
+
+        /// Viscosity of stock components at surface conditions.
+        /// \return Array of 2 viscosity values.
+        virtual const double* viscosity() const = 0;
+
+
+
+        // ------ Relative permeability ------
+
+        /// Relative permeabilities for all phases.
+        /// \param[in]  sw     Array of n water saturation values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the saturation values.
+        /// \return            An std::vector with 2 elements, each an array of n relperm values,
+        ///                    containing krw, kro. Use PhaseIndex for indexing into the result.
+        virtual
+        std::vector<V> relperm(const V& sw,
+                               const V& so,
+                               const Cells& cells) const = 0;
+
+        /// Relative permeabilities for all phases.
+        /// \param[in]  sw     Array of n water saturation values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the saturation values.
+        /// \return            An std::vector with 2 elements, each an array of n relperm values,
+        ///                    containing krw, kro. Use PhaseIndex for indexing into the result.
+        virtual
+        std::vector<ADB> relperm(const ADB& sw,
+                                 const ADB& so,
+                                 const Cells& cells) const = 0;
+
+    };
+} // namespace Opm
+
+#endif// OPM_INCOMPPROPSADINTERFACE_HEADER_INCLUDED

opm/polymer/fullyimplicit/PolymerPropsAd.cpp

@@ -0,0 +1,216 @@
+
+#include <cmath>
+#include <vector>
+#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
+#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
+#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+
+namespace Opm {
+
+
+    typedef PolymerPropsAd::ADB ADB;
+    typedef PolymerPropsAd::V V;
+
+/*
+    PolymerPropsAd::PolymerPropsAd()
+    {
+    }
+
+
+    PolymerPropsAd::PolymerPropsAd(const int num_cells,
+                                   const double c_max,
+                                   const double mix_param,
+                                   const double std::vector<double>& c_vals_visc,
+                                   const double std::vector<double>& visc_mult_vals)
+        : num_cells_ (num_cells)
+        , c_max_ (c_max)
+        , mix_param_(mix_param)
+        , c_vals_visc_ (c_vals_visc)
+        , visc_mult_vals_ (visc_mult_vals)
+    {
+    }
+    
+
+    double PolymerProsAd::num_cells() const
+    {
+        return num__cells_;
+    }
+
+
+
+
+    double PolymerPropsAd::cMax() const
+    {
+        return c_max_;
+    }
+
+
+
+    double PolymerPropsAd::mixParam() const
+    {
+        return mix_param_;
+    }
+
+
+
+    V PolymerPropsAd::muM(const V& c,
+                          const double* visc) const
+    {
+        const int nc = num_cells();
+        assert(nc == c.size());
+        std::vector<double> mu(nc);
+
+        for (int i = 0; i < nc; ++i) {
+            mu[i] = Opm::linearInterpolation(c_vals_visc_, visc_mult_vals_, c(i));
+        }
+
+        const V muM = Eigen::Map<const V>(&mu[0], nc);
+        
+        const double mu_w = visc[0];
+
+        return muM * mu_w;
+    }
+
+   
+
+    ADB PolymerPropsAd::muM(const ADB& c,
+                            const double* visc) const
+    {
+        const int nc = num_cells();
+        assert(nc == c.size());
+        
+        V mu_m = muM(c.value());
+
+        std::vector<double> dmu_dc(nc);
+
+        for (int i = 0; i < nc; ++i) {
+            dmu_dc[i] = Opm::linearInterpolationDerivative(c_vals_visc_, visc_mult_vals_, c.value()(i));
+        }      
+        
+        const V dmu = Eigen::Map<const V>(&dmu_dc[0], nc);
+
+        ADB::M dmu_diag = spdiag(dmu);
+        const int num_blocks = c.numBlocks();
+        std::vector<ADB::M> jacs(num_blocks);
+
+        for (int block = 0; block < num_blocks; ++block) {
+            jacs[block] = dmu_diag * c.derivative()[block];
+        }
+        
+        const double mu_w = visc[0]
+
+        return ADB::function(mu_m, jacs) * mu_w;
+    }
+
+
+
+    V PolymerPropsAd::ToddLongstaff(const double mix_param,
+                                    const V& muM,
+                                    const V& mu) const
+    {
+        const int nc = num_cells();
+        assert(nc == muM.size());
+
+        std::vector<double> mueff(nc);
+        const double omega = mix_param;
+
+        for (int i = 0; i < nc; ++i) {
+            mueff[i] = std::pow(muM(i),omega) * std::pow(mu(i), 1. - omega);
+        }
+
+        const V muEff = Eigen::Map<const V>(&mueff[0], nc);
+
+        return muEff;
+    }
+
+
+    ADB PolymerPropsAd::ToddLongstaff(const double mix_param,
+                                      const ADB& muM,
+                                      const ADB& mu) const
+    {
+        const int nc = num_cells();
+
+    }
+
+
+    V PolymerPropsAd::muPolyEff(const double mix_param,
+                                const V& muM,
+                                const V& muPoly) const
+    {
+    
+        return ToddLongstaff(mix_param, muM, muPoly);
+    }
+
+    V PolymerPropsAd::muWatEff(const double mix_param,
+                               const std::vecotr<double>& c_max,
+                               const V& c,
+                               const V& muM,
+                               const V& muWat,
+                               const V& muPolyEff) const
+    {
+        const int nc = num_cells();
+        assert(nc == c.size());
+        V muWate = ToddLongstaff(mix_param, muM, muWat);
+
+//        V cmax = V::Constant(nc, 1, c_max);
+        const V cmax = Eigen::Map<const V>(&c_max[0], nc);
+        const V one = V::Ones(nc, 1);
+        V inv_muWatEff = (one - c / cmax) / muWate + c / cmax / muPolyEff;
+
+        V muWatEff = one / inv_muWatEff;
+
+        return muWatEff;
+
+    }
+*/
+    V PolymerPropsAd::effectiveInvWaterVisc(const V& c,
+                                            const double* visc) const
+    {
+        const int nc = c.size();
+        V inv_mu_w_eff(n);
+        for (int i = 0; i < nc; ++i) {
+            double im;
+            polymer_props_.effectiveInvVisc(c(i), visc, im);
+            inv_mu_w_eff(i) = im;
+        }
+
+        return inv_mu_w_eff;
+    }
+
+
+
+
+
+
+    ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c,
+	                    				      const double* visc)
+    {
+	    const int n = c.size();
+    	V inv_mu_w_eff(n);
+    	V dinv_mu_w_eff(n);
+    	for (int i = 0; i < n; ++i) {
+    	    double im, dim;
+    	    polymer_props_.effectiveInvViscWithDer(c.value()(i), visc, im, dim);
+    	    inv_mu_w_eff(i) = im;
+    	    dinv_mu_w_eff(i) = dim;
+    	}
+            ADB::M dim_diag = spdiag(dinv_mu_w_eff);
+            const int num_blocks = c.numBlocks();
+            std::vector<ADB::M> jacs(num_blocks);
+            for (int block = 0; block < num_blocks; ++block) {
+                jacs[block] = dim_diag * c.derivative()[block];
+            }
+        return ADB::function(inv_mu_w_eff, jacs);
+    }
+
+
+
+
+
+    V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const
+    {
+        const int nc
+       
+    }
+
+} // namespace Opm

opm/polymer/fullyimplicit/PolymerPropsAd.hpp

@@ -0,0 +1,89 @@
+#ifndef OPM_POLYMERPROPSAD_HEADED_INLCUDED
+#define OPM_POLYMERPROPSAD_HEADED_INLCUDED
+
+
+#include <cmath>
+#include <vector>
+#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
+#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
+#include <opm/polymer/polymerProperties.hpp>
+namespace Opm {
+
+    
+    class PolymerPropsAd 
+    {
+    public:
+/*        PolymerPropsAd(const int num_cells;
+                       const double mix_param,
+                       const std::vector<double>& c_max,
+                       const std::vector<double>& c_vals_visc,
+                       const std::vector<double>& visc_mult_vals);
+        
+
+        double num_cells() const;
+        double cMax() const;
+        double mixParam() const;
+
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+
+
+        V muM(const V& c,
+              const double* visc) const;
+        ADB muM(const ADB& c
+                const double* visc) const;
+
+        V ToddLongstaff(const double mix_param,
+                        const V& muM,
+                        const V& mu) const;
+
+        ADB ToddLongstaff(const double mix_param,
+                          const ADB& muM,
+                          const ADB& mu) const;
+       
+        V muPolyEff(const double mix_param,
+                    const V& muM,
+                    const V& muPoly) const;
+    
+        ADB muPolyEff(const double mix_param,
+                      const ADB& muM,
+                      const ADB& muPoly) const;
+
+        V muWatEff(const double mix_param,
+                   const std::vector<double>& c_max,
+                   const V& c,
+                   const V& muM,
+                   const V& muWat,
+                   const V& muPolyEff) const;
+
+        ADB muWatEff(const double mix_param,
+                     const std::vector<double>& c_max,
+                     const ADB& c,
+                     const ADB& muM,
+                     const ADB& muWat,
+                     const ADB& muPolyEff) const;
+*/
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        PolymerPropsAd(const PolymerPropperties& polyprops);
+        ~PolymerPropsAd();
+    private:
+        const PolymerProperties polyprops_;
+    };
+
+    
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+#endif// OPM_POLYMERPROPSAD_HEADED_INLCUDED

opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.cpp

@@ -0,0 +1,261 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/autodiff/polymer/SimulatorFullyImplicitTwophase.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+#include <opm/autodiff/polymer/FullyImplicitTwoPhaseSolver.hpp>
+#include <opm/autodiff/polymer/IncompPropsAdInterface.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/wells.h>
+#include <opm/core/pressure/flow_bc.h>
+
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/io/vtk/writeVtkData.hpp>
+
+
+#include <opm/core/utility/miscUtilities.hpp>
+
+#include <opm/core/grid/ColumnExtract.hpp>
+#include <opm/core/simulator/TwophaseState.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/scoped_ptr.hpp>
+#include <boost/lexical_cast.hpp>
+
+#include <numeric>
+#include <fstream>
+#include <iostream>
+#include <Eigen/Eigen>
+namespace Opm
+{
+
+    class SimulatorFullyImplicitTwophase::Impl
+    {
+    public:
+        Impl(const parameter::ParameterGroup& param,
+             const UnstructuredGrid& grid,
+             const IncompPropsAdInterface& props,
+             LinearSolverInterface& linsolver,
+             std::vector<double>& src);
+
+        SimulatorReport run(SimulatorTimer& timer,
+                            TwophaseState& state,
+                            std::vector<double>& src);
+
+    private:
+
+        // Parameters for output.
+        bool output_;
+        bool output_vtk_;
+        std::string output_dir_;
+        int output_interval_;
+        // Parameters for well control
+        // Observed objects.
+        const UnstructuredGrid& grid_;
+        const IncompPropsAdInterface& props_;
+        // Solvers
+        FullyImplicitTwoPhaseSolver solver_;
+        // Misc. data
+        std::vector<int> allcells_;
+    };
+
+
+
+
+    SimulatorFullyImplicitTwophase::SimulatorFullyImplicitTwophase(const parameter::ParameterGroup& param,
+                                                                   const UnstructuredGrid& grid,
+                                                                   const IncompPropsAdInterface& props,
+                                                                   LinearSolverInterface& linsolver,
+                                                                   std::vector<double>& src)
+    {
+        pimpl_.reset(new Impl(param, grid, props, linsolver, src));
+    }
+
+
+
+
+
+    SimulatorReport SimulatorFullyImplicitTwophase::run(SimulatorTimer& timer,
+                                                        TwophaseState& state,
+                                                        std::vector<double>& src)
+    {
+        return pimpl_->run(timer, state, src);
+    }
+
+
+
+    static void outputStateVtk(const UnstructuredGrid& grid,
+                               const Opm::TwophaseState& state,
+                               const int step,
+                               const std::string& output_dir)
+    {
+        // Write data in VTK format.
+        std::ostringstream vtkfilename;
+        vtkfilename << output_dir << "/vtk_files";
+        boost::filesystem::path fpath(vtkfilename.str());
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
+        std::ofstream vtkfile(vtkfilename.str().c_str());
+        if (!vtkfile) {
+            OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str());
+        }
+        Opm::DataMap dm;
+        dm["saturation"] = &state.saturation();
+        dm["pressure"] = &state.pressure();
+        std::vector<double> cell_velocity;
+        Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
+        dm["velocity"] = &cell_velocity;
+        Opm::writeVtkData(grid, dm, vtkfile);
+    }
+/*
+    static void outputWaterCut(const Opm::Watercut& watercut,
+                               const std::string& output_dir)
+    {
+        // Write water cut curve.
+        std::string fname = output_dir  + "/watercut.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        watercut.write(os);
+    }
+    
+   */ 
+    
+    
+    
+    
+    SimulatorFullyImplicitTwophase::Impl::Impl(const parameter::ParameterGroup& param,
+                                               const UnstructuredGrid& grid,
+                                               const IncompPropsAdInterface& props,
+                                               LinearSolverInterface& linsolver,
+                                               std::vector<double>& src)
+        : grid_(grid),
+          props_(props),
+          solver_(grid_, props_, linsolver)
+    {
+        // For output.
+        output_ = param.getDefault("output", true);
+        if (output_) {
+            output_vtk_ = param.getDefault("output_vtk", true);
+            output_dir_ = param.getDefault("output_dir", std::string("output"));
+            // Ensure that output dir exists
+            boost::filesystem::path fpath(output_dir_);
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+            }
+            output_interval_ = param.getDefault("output_interval", 1);
+        }
+
+        // Misc init.
+        const int num_cells = grid.number_of_cells;
+        allcells_.resize(num_cells);
+        for (int cell = 0; cell < num_cells; ++cell) {
+            allcells_[cell] = cell;
+        }
+    }
+
+
+
+
+    SimulatorReport SimulatorFullyImplicitTwophase::Impl::run(SimulatorTimer& timer,
+                                                              TwophaseState& state,
+                                                              std::vector<double>& src)
+    {
+        // Initialisation.
+        std::vector<double> porevol;
+        computePorevolume(grid_, props_.porosity(), porevol);
+        // const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
+        std::vector<double> initial_porevol = porevol;
+
+        // Main simulation loop.
+        Opm::time::StopWatch solver_timer;
+        double stime = 0.0;
+        Opm::time::StopWatch step_timer;
+        Opm::time::StopWatch total_timer;
+        total_timer.start();
+        // These must be changed for three-phase.
+        std::vector<double> fractional_flows;
+        std::fstream tstep_os;
+        if (output_) {
+            std::string filename = output_dir_ + "/step_timing.param";
+            tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
+        }
+        for (; !timer.done(); ++timer) {
+            // Report timestep and (optionally) write state to disk.
+            step_timer.start();
+            timer.report(std::cout);
+            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
+                if (output_vtk_) {
+                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+                }
+            }
+   
+            SimulatorReport sreport;
+
+            // Run solver.
+            solver_timer.start();
+            std::vector<double> initial_pressure = state.pressure();
+            solver_.step(timer.currentStepLength(), state, src);
+
+            // Stop timer and report.
+            solver_timer.stop();
+            const double st = solver_timer.secsSinceStart();
+            std::cout << "Fully implicit solver took:  " << st << " seconds." << std::endl;
+            stime += st;
+            sreport.pressure_time = st;
+
+            sreport.total_time =  step_timer.secsSinceStart();
+            if (output_) {
+                sreport.reportParam(tstep_os);
+            }
+        }
+
+        if (output_) {
+            if (output_vtk_) {
+                outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+            }
+     //       outputWaterCut(watercut, output_dir_);
+            tstep_os.close();
+        }
+
+        total_timer.stop();
+
+        SimulatorReport report;
+        report.pressure_time = stime;
+        report.transport_time = 0.0;
+        report.total_time = total_timer.secsSinceStart();
+        return report;
+    }
+
+
+} // namespace Opm

opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.hpp

@@ -0,0 +1,85 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_SIMULATORFULLYIMPLICITTWOPHASE_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITTWOPHASE_HEADER_INCLUDED
+
+#include <boost/shared_ptr.hpp>
+#include <vector>
+
+struct UnstructuredGrid;
+
+namespace Opm
+{
+    namespace parameter { class ParameterGroup; }
+    class IncompPropsAdInterface;
+    class LinearSolverInterface;
+    class SimulatorTimer;
+    class TwophaseState;
+    struct SimulatorReport;
+
+    /// Class collecting all necessary components for a two-phase simulation.
+    class SimulatorFullyImplicitTwophase
+    {
+    public:
+        /// Initialise from parameters and objects to observe.
+        /// \param[in] param       parameters, this class accepts the following:
+        ///     parameter (default)            effect
+        ///     -----------------------------------------------------------
+        ///     output (true)                  write output to files?
+        ///     output_dir ("output")          output directoty
+        ///     output_interval (1)            output every nth step
+        ///     nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
+        ///     nl_pressure_change_tolerance (1.0)   pressure solver change tolerance (in Pascal)
+        ///     nl_pressure_maxiter (10)       max nonlinear iterations in pressure
+        ///     nl_maxiter (30)                max nonlinear iterations in transport
+        ///     nl_tolerance (1e-9)            transport solver absolute residual tolerance
+        ///     num_transport_substeps (1)     number of transport steps per pressure step
+        ///     use_segregation_split (false)  solve for gravity segregation (if false,
+        ///                                    segregation is ignored).
+        ///
+        /// \param[in] grid          grid data structure
+        /// \param[in] props         fluid and rock properties
+        /// \param[in] linsolver     linear solver
+        SimulatorFullyImplicitTwophase(const parameter::ParameterGroup& param,
+                                       const UnstructuredGrid& grid,
+                                       const IncompPropsAdInterface& props,
+                                       LinearSolverInterface& linsolver,
+                                       std::vector<double>& src);
+
+        /// Run the simulation.
+        /// This will run succesive timesteps until timer.done() is true. It will
+        /// modify the reservoir and well states.
+        /// \param[in,out] timer       governs the requested reporting timesteps
+        /// \param[in,out] state       state of reservoir: pressure, fluxes
+        /// \param[in,out] well_state  state of wells: bhp, perforation rates
+        /// \return                    simulation report, with timing data
+        SimulatorReport run(SimulatorTimer& timer,
+                            TwophaseState& state,
+                            std::vector<double>& src);
+
+    private:
+        class Impl;
+        // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl.
+        boost::shared_ptr<Impl> pimpl_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED