Add an example program of FIBOS that uses CpGrid.

opm/autodiff/AutoDiffHelpers.hpp

@@ -59,7 +59,8 @@ struct HelperOps
     M fulldiv;
 
     /// Constructs all helper vectors and matrices.
-    HelperOps(const UnstructuredGrid& grid)
+    template<class Grid>
+    HelperOps(const Grid& grid)
     {
         using namespace AutoDiffGrid;
         const int nc = numCells(grid);
@@ -93,7 +94,7 @@ struct HelperOps
         div = ngrad.transpose();
         std::vector<Tri> fullngrad_tri;
         fullngrad_tri.reserve(2*nf);
-        ADFaceCellTraits<UnstructuredGrid>::Type nb=faceCells(grid);
+        typename ADFaceCellTraits<Grid>::Type nb=faceCells(grid);
         for (int i = 0; i < nf; ++i) {
             if (nb(i,0) >= 0) {
                 fullngrad_tri.emplace_back(i, nb(i,0), 1.0);
@@ -118,14 +119,15 @@ struct HelperOps
     public:
         typedef AutoDiffBlock<Scalar> ADB;
 
-        UpwindSelector(const UnstructuredGrid& g,
+        template<class Grid>
+        UpwindSelector(const Grid& g,
                        const HelperOps&        h,
                        const typename ADB::V&  ifaceflux)
         {
             using namespace AutoDiffGrid;
             typedef HelperOps::IFaces::Index IFIndex;
             const IFIndex nif = h.internal_faces.size();
-            ADFaceCellTraits<UnstructuredGrid>::Type
+            typename ADFaceCellTraits<Grid>::Type
                 face_cells = faceCells(g);
             assert(nif == ifaceflux.size());
 

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -45,9 +45,12 @@ namespace Opm {
     ///
     /// It uses automatic differentiation via the class AutoDiffBlock
     /// to simplify assembly of the jacobian matrix.
+    template<class T>
     class FullyImplicitBlackoilSolver
     {
     public:
+        /// \brief The type of the grid that we use.
+        typedef T Grid;
         /// Construct a solver. It will retain references to the
         /// arguments of this functions, and they are expected to
         /// remain in scope for the lifetime of the solver.
@@ -57,7 +60,7 @@ namespace Opm {
         /// \param[in] rock_comp_props  if non-null, rock compressibility properties
         /// \param[in] wells            well structure
         /// \param[in] linsolver        linear solver
-        FullyImplicitBlackoilSolver(const UnstructuredGrid&         grid ,
+        FullyImplicitBlackoilSolver(const Grid&                     grid ,
                                     const BlackoilPropsAdInterface& fluid,
                                     const DerivedGeology&           geo  ,
                                     const RockCompressibility*      rock_comp_props,
@@ -118,7 +121,7 @@ namespace Opm {
                Gas   = BlackoilPropsAdInterface::Gas  };
 
         // Member data
-        const UnstructuredGrid&         grid_;
+        const Grid&         grid_;
         const BlackoilPropsAdInterface& fluid_;
         const DerivedGeology&           geo_;
         const RockCompressibility*      rock_comp_props_;
@@ -250,5 +253,6 @@ namespace Opm {
     };
 } // namespace Opm
 
+#include "FullyImplicitBlackoilSolver_impl.hpp"
 
 #endif // OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -74,22 +74,23 @@ namespace {
         return all_cells;
     }
 
-    template <class GeoProps>
+    template <class GeoProps, class Grid>
     AutoDiffBlock<double>::M
-    gravityOperator(const UnstructuredGrid& grid,
+    gravityOperator(const Grid&             grid,
                     const HelperOps&        ops ,
                     const GeoProps&         geo )
     {
         using namespace Opm::AutoDiffGrid;
         const int nc = numCells(grid);
-        SparseTableView c2f = cell2Faces(grid);
+        typedef typename Opm::UgGridHelpers::Cell2FacesTraits<Grid>::Type Cell2Faces;
+        Cell2Faces c2f = cell2Faces(grid);
 
         std::vector<int> f2hf(2 * numFaces(grid), -1);
-        Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
+        typename ADFaceCellTraits<Grid>::Type
             face_cells = faceCells(grid);
         for (int c = 0, i = 0; c < nc; ++c) {
-            typename SparseTableView::row_type faces=c2f[c];
-            typedef typename SparseTableView::row_type::iterator Iter;
+            typename Cell2Faces::row_type faces=c2f[c];
+            typedef typename Cell2Faces::row_type::iterator Iter;
             for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f) {
                 const int p = 0 + (face_cells(*f, 0) != c);
 
@@ -128,8 +129,8 @@ namespace {
     }
 
 
-
-    V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav)
+    template<class Grid>
+    V computePerfPress(const Grid& grid, const Wells& wells, const V& rho, const double grav)
     {
         using namespace Opm::AutoDiffGrid;
         const int nw = wells.number_of_wells;
@@ -194,9 +195,9 @@ namespace {
 
 
 
-
-    FullyImplicitBlackoilSolver::
-    FullyImplicitBlackoilSolver(const UnstructuredGrid&         grid ,
+    template<class T>
+    FullyImplicitBlackoilSolver<T>::
+    FullyImplicitBlackoilSolver(const Grid&                     grid ,
                                 const BlackoilPropsAdInterface& fluid,
                                 const DerivedGeology&           geo  ,
                                 const RockCompressibility*      rock_comp_props,
@@ -215,7 +216,7 @@ namespace {
         , wops_  (wells)
         , grav_  (gravityOperator(grid_, ops_, geo_))
         , rq_    (fluid.numPhases())
-        , phaseCondition_(grid.number_of_cells)
+        , phaseCondition_(AutoDiffGrid::numCells(grid))
         , residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
                         ADB::null(),
                         ADB::null() } )
@@ -225,9 +226,9 @@ namespace {
 
 
 
-
+    template<class T>
     void
-    FullyImplicitBlackoilSolver::
+    FullyImplicitBlackoilSolver<T>::
     step(const double   dt,
          BlackoilState& x ,
          WellState&     xw)
@@ -278,7 +279,8 @@ namespace {
 
 
 
-    FullyImplicitBlackoilSolver::ReservoirResidualQuant::ReservoirResidualQuant()
+    template<class T>
+    FullyImplicitBlackoilSolver<T>::ReservoirResidualQuant::ReservoirResidualQuant()
         : accum(2, ADB::null())
         , mflux(   ADB::null())
         , b    (   ADB::null())
@@ -291,7 +293,8 @@ namespace {
 
 
 
-    FullyImplicitBlackoilSolver::SolutionState::SolutionState(const int np)
+    template<class T>
+    FullyImplicitBlackoilSolver<T>::SolutionState::SolutionState(const int np)
         : pressure  (    ADB::null())
         , saturation(np, ADB::null())
         , rs        (    ADB::null())
@@ -305,7 +308,8 @@ namespace {
 
 
 
-    FullyImplicitBlackoilSolver::
+    template<class T>
+    FullyImplicitBlackoilSolver<T>::
     WellOps::WellOps(const Wells& wells)
         : w2p(wells.well_connpos[ wells.number_of_wells ],
               wells.number_of_wells)
@@ -336,8 +340,9 @@ namespace {
 
 
 
-    FullyImplicitBlackoilSolver::SolutionState
-    FullyImplicitBlackoilSolver::constantState(const BlackoilState& x,
+    template<class T>
+    typename FullyImplicitBlackoilSolver<T>::SolutionState
+    FullyImplicitBlackoilSolver<T>::constantState(const BlackoilState& x,
                                                const WellState&     xw)
     {
         using namespace Opm::AutoDiffGrid;
@@ -433,8 +438,9 @@ namespace {
 
 
 
-    FullyImplicitBlackoilSolver::SolutionState
-    FullyImplicitBlackoilSolver::variableState(const BlackoilState& x,
+    template<class T>
+    typename FullyImplicitBlackoilSolver<T>::SolutionState
+    FullyImplicitBlackoilSolver<T>::variableState(const BlackoilState& x,
                                                const WellState&     xw)
     {
         using namespace Opm::AutoDiffGrid;
@@ -580,8 +586,9 @@ namespace {
 
 
 
+    template<class T>
     void
-    FullyImplicitBlackoilSolver::computeAccum(const SolutionState& state,
+    FullyImplicitBlackoilSolver<T>::computeAccum(const SolutionState& state,
                                               const int            aix  )
     {
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
@@ -621,8 +628,9 @@ namespace {
 
 
 
+    template<class T>
     void
-    FullyImplicitBlackoilSolver::
+    FullyImplicitBlackoilSolver<T>::
     assemble(const V&             pvdt,
              const BlackoilState& x   ,
              const WellState&     xw  )
@@ -834,7 +842,8 @@ namespace {
 
 
 
-    V FullyImplicitBlackoilSolver::solveJacobianSystem() const
+    template<class T>
+    V FullyImplicitBlackoilSolver<T>::solveJacobianSystem() const
     {
         const int np = fluid_.numPhases();
         ADB mass_res = residual_.mass_balance[0];
@@ -886,7 +895,8 @@ namespace {
 
 
 
-    void FullyImplicitBlackoilSolver::updateState(const V& dx,
+    template<class T>
+    void FullyImplicitBlackoilSolver<T>::updateState(const V& dx,
                                                   BlackoilState& state,
                                                   WellState& well_state)
     {
@@ -1136,8 +1146,9 @@ namespace {
 
 
 
+    template<class T>
     std::vector<ADB>
-    FullyImplicitBlackoilSolver::computeRelPerm(const SolutionState& state) const
+    FullyImplicitBlackoilSolver<T>::computeRelPerm(const SolutionState& state) const
     {
         using namespace Opm::AutoDiffGrid;
         const int               nc   = numCells(grid_);
@@ -1162,8 +1173,9 @@ namespace {
     }
 
 
+    template<class T>
     std::vector<ADB>
-    FullyImplicitBlackoilSolver::computePressures(const SolutionState& state) const
+    FullyImplicitBlackoilSolver<T>::computePressures(const SolutionState& state) const
     {
         using namespace Opm::AutoDiffGrid;
         const int               nc   = numCells(grid_);
@@ -1204,8 +1216,9 @@ namespace {
 
 
 
+    template<class T>
     std::vector<ADB>
-    FullyImplicitBlackoilSolver::computeRelPermWells(const SolutionState& state,
+    FullyImplicitBlackoilSolver<T>::computeRelPermWells(const SolutionState& state,
                                                      const DataBlock& well_s,
                                                      const std::vector<int>& well_cells) const
     {
@@ -1235,8 +1248,9 @@ namespace {
 
 
 
+    template<class T>
     void
-    FullyImplicitBlackoilSolver::computeMassFlux(const int               actph ,
+    FullyImplicitBlackoilSolver<T>::computeMassFlux(const int               actph ,
                                                  const V&                transi,
                                                  const ADB&              kr    ,
                                                  const ADB&              phasePressure,
@@ -1276,8 +1290,9 @@ namespace {
 
 
 
+    template<class T>
     double
-    FullyImplicitBlackoilSolver::residualNorm() const
+    FullyImplicitBlackoilSolver<T>::residualNorm() const
     {
         double r = 0;
         for (std::vector<ADB>::const_iterator
@@ -1297,8 +1312,9 @@ namespace {
 
 
 
+    template<class T>
     ADB
-    FullyImplicitBlackoilSolver::fluidViscosity(const int               phase,
+    FullyImplicitBlackoilSolver<T>::fluidViscosity(const int               phase,
                                                 const ADB&              p    ,
                                                 const ADB&              rs   ,
                                                 const ADB&              rv   ,
@@ -1322,8 +1338,9 @@ namespace {
 
 
 
+    template<class T>
     ADB
-    FullyImplicitBlackoilSolver::fluidReciprocFVF(const int               phase,
+    FullyImplicitBlackoilSolver<T>::fluidReciprocFVF(const int               phase,
                                                   const ADB&              p    ,
                                                   const ADB&              rs   ,
                                                   const ADB&              rv   ,
@@ -1347,8 +1364,9 @@ namespace {
 
 
 
+    template<class T>
     ADB
-    FullyImplicitBlackoilSolver::fluidDensity(const int               phase,
+    FullyImplicitBlackoilSolver<T>::fluidDensity(const int               phase,
                                               const ADB&              p    ,
                                               const ADB&              rs   ,
                                               const ADB&              rv   ,
@@ -1373,8 +1391,9 @@ namespace {
 
 
 
+    template<class T>
     V
-    FullyImplicitBlackoilSolver::fluidRsSat(const V&                p,
+    FullyImplicitBlackoilSolver<T>::fluidRsSat(const V&                p,
                                             const std::vector<int>& cells) const
     {
         return fluid_.rsSat(p, cells);
@@ -1384,15 +1403,17 @@ namespace {
 
 
 
+    template<class T>
     ADB
-    FullyImplicitBlackoilSolver::fluidRsSat(const ADB&              p,
+    FullyImplicitBlackoilSolver<T>::fluidRsSat(const ADB&              p,
                                             const std::vector<int>& cells) const
     {
         return fluid_.rsSat(p, cells);
     }
 
+    template<class T>
     V
-    FullyImplicitBlackoilSolver::fluidRvSat(const V&                p,
+    FullyImplicitBlackoilSolver<T>::fluidRvSat(const V&                p,
                                             const std::vector<int>& cells) const
     {
         return fluid_.rvSat(p, cells);
@@ -1402,8 +1423,9 @@ namespace {
 
 
 
+    template<class T>
     ADB
-    FullyImplicitBlackoilSolver::fluidRvSat(const ADB&              p,
+    FullyImplicitBlackoilSolver<T>::fluidRvSat(const ADB&              p,
                                             const std::vector<int>& cells) const
     {
         return fluid_.rvSat(p, cells);
@@ -1411,8 +1433,9 @@ namespace {
 
 
 
+    template<class T>
     ADB
-    FullyImplicitBlackoilSolver::poroMult(const ADB& p) const
+    FullyImplicitBlackoilSolver<T>::poroMult(const ADB& p) const
     {
         const int n = p.size();
         if (rock_comp_props_ && rock_comp_props_->isActive()) {
@@ -1438,8 +1461,9 @@ namespace {
 
 
 
+    template<class T>
     ADB
-    FullyImplicitBlackoilSolver::transMult(const ADB& p) const
+    FullyImplicitBlackoilSolver<T>::transMult(const ADB& p) const
     {
         const int n = p.size();
         if (rock_comp_props_ && rock_comp_props_->isActive()) {
@@ -1463,8 +1487,9 @@ namespace {
 
 
     /*
+    template<class T>
     void
-    FullyImplicitBlackoilSolver::
+    FullyImplicitBlackoilSolver<T>::
     classifyCondition(const SolutionState&        state,
                       std::vector<PhasePresence>& cond ) const
     {
@@ -1504,10 +1529,12 @@ namespace {
     } */
 
 
+    template<class T>
     void
-    FullyImplicitBlackoilSolver::classifyCondition(const BlackoilState& state)
+    FullyImplicitBlackoilSolver<T>::classifyCondition(const BlackoilState& state)
     {
-        const int nc = grid_.number_of_cells;
+        using namespace Opm::AutoDiffGrid;
+        const int nc = numCells(grid_);
         const int np = state.numPhases();
 
         const PhaseUsage& pu = fluid_.phaseUsage();

opm/autodiff/GeoProps.hpp

@@ -21,8 +21,9 @@
 #define OPM_GEOPROPS_HEADER_INCLUDED
 
 #include <opm/core/grid.h>
-#include <opm/core/pressure/tpfa/trans_tpfa.h>
 #include <opm/autodiff/GridHelpers.hpp>
+//#include <opm/core/pressure/tpfa/trans_tpfa.h>
+#include <opm/core/pressure/tpfa/TransTpfa.hpp>
 #include <Eigen/Eigen>
 
 namespace Opm
@@ -40,8 +41,8 @@ namespace Opm
 
         /// Construct contained derived geological properties
         /// from grid and property information.
-        template <class Props>
-        DerivedGeology(const UnstructuredGrid& grid,
+        template <class Props, class Grid>
+        DerivedGeology(const Grid&             grid,
                        const Props&            props ,
                        const double*           grav = 0)
             : pvol_ (Opm::AutoDiffGrid::numCells(grid))
@@ -57,8 +58,8 @@ namespace Opm
                            std::multiplies<double>());
 
             // Transmissibility
-            Vector htrans(grid.cell_facepos[nc]);
-            UnstructuredGrid* ug = const_cast<UnstructuredGrid*>(& grid);
+            Vector htrans(numCellFaces(grid));
+            Grid* ug = const_cast<Grid*>(& grid);
             tpfa_htrans_compute(ug, props.permeability(), htrans.data());
             tpfa_trans_compute (ug, htrans.data()     , trans_.data());
 
@@ -72,13 +73,14 @@ namespace Opm
             std::fill(gravity_, gravity_ + 3, 0.0);
             if (grav != 0) {
                 const typename Vector::Index nd = dimensions(grid);
-                SparseTableView c2f=cell2Faces(grid);
+                typedef typename ADCell2FacesTraits<Grid>::Type Cell2Faces;
+                Cell2Faces c2f=cell2Faces(grid);
 
                 for (typename Vector::Index c = 0; c < nc; ++c) {
                     const double* const cc = cellCentroid(grid, c);
 
-                    typename SparseTableView::row_type faces=c2f[c];
-                    typedef SparseTableView::row_type::iterator Iter;
+                    typename Cell2Faces::row_type faces=c2f[c];
+                    typedef typename Cell2Faces::row_type::iterator Iter;
 
                     for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f) {
                         const double* const fc = faceCentroid(grid, *f);

opm/autodiff/GridHelpers.cpp

@@ -17,7 +17,9 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
+
 #include "config.h"
+
 #include <opm/autodiff/GridHelpers.hpp>
 namespace Opm
 {
@@ -112,10 +114,14 @@ void extractInternalFaces(const UnstructuredGrid& grid,
         }
     }
 }
+} // end namespace AutoDiffHelpers
 
 #ifdef HAVE_DUNE_CORNERPOINT
 // Interface functions using CpGrid
 
+namespace UgGridHelpers
+{
+
 int numCells(const Dune::CpGrid& grid)
 {
     return grid.numCells();
@@ -131,6 +137,11 @@ int dimensions(const Dune::CpGrid&)
     return Dune::CpGrid::dimension;
 }
 
+int numCellFaces(const Dune::CpGrid& grid)
+{
+    return grid.numCellFaces();    
+}
+
 const int* cartDims(const Dune::CpGrid& grid)
 {
     return &(grid.logicalCartesianSize()[0]);
@@ -141,11 +152,51 @@ const int*  globalCell(const Dune::CpGrid& grid)
     return &(grid.globalCell()[0]);
 }
 
-FaceCellsContainerProxy faceCells(const Dune::CpGrid& grid)
+CellCentroidTraits<Dune::CpGrid>::IteratorType
+beginCellCentroids(const Dune::CpGrid& grid)
 {
-    return FaceCellsContainerProxy(&grid);
+    return CellCentroidTraits<Dune::CpGrid>::IteratorType(grid, 0);
 }
 
+double cellCentroidCoordinate(const Dune::CpGrid& grid, int cell_index,
+                              int coordinate)
+{
+    return grid.cellCentroid(cell_index)[coordinate];
+}
+
+FaceCentroidTraits<Dune::CpGrid>::IteratorType
+beginFaceCentroids(const Dune::CpGrid& grid)
+{
+    return FaceCentroidTraits<Dune::CpGrid>::IteratorType(grid, 0);
+}
+
+FaceCentroidTraits<Dune::CpGrid>::ValueType
+faceCentroid(const Dune::CpGrid& grid, int face_index)
+{
+    return grid.faceCentroid(face_index);
+}
+
+Opm::AutoDiffGrid::Cell2FacesContainer cell2Faces(const Dune::CpGrid& grid)
+{
+    return Opm::AutoDiffGrid::Cell2FacesContainer(&grid);
+}
+
+FaceCellTraits<Dune::CpGrid>::Type
+faceCells(const Dune::CpGrid& grid)
+{
+    return Opm::AutoDiffGrid::FaceCellsContainerProxy(&grid);
+}
+
+const double* faceNormal(const Dune::CpGrid& grid, int face_index)
+{
+    return &(grid.faceNormal(face_index)[0]);
+}
+
+} // end namespace UgGridHelpers
+
+namespace AutoDiffGrid
+{
+
 
 Eigen::Array<double, Eigen::Dynamic, 1>
 cellCentroidsZ(const Dune::CpGrid& grid)
@@ -169,11 +220,6 @@ const double* faceCentroid(const Dune::CpGrid& grid, int face_index)
     return &(grid.faceCentroid(face_index)[0]);
 }
 
-Cell2FacesContainer cell2Faces(const Dune::CpGrid& grid)
-{
-    return Cell2FacesContainer(&grid);
-}
-
 double cellVolume(const  Dune::CpGrid& grid, int cell_index)
 {
     return grid.cellVolume(cell_index);

opm/autodiff/GridHelpers.hpp

@@ -20,6 +20,8 @@
 #ifndef OPM_GRIDHELPERS_HEADER_INCLUDED
 #define OPM_GRIDHELPERS_HEADER_INCLUDED
 
+#include <functional>
+
 #include <boost/range/iterator_range.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/grid/GridHelpers.hpp>
@@ -33,18 +35,12 @@
 namespace Opm
 {
 
-namespace AutoDiffGrid
+namespace UgGridHelpers
 {
 
-using Opm::UgGridHelpers::SparseTableView;
-using Opm::UgGridHelpers::numCells;
-using Opm::UgGridHelpers::numFaces;
-using Opm::UgGridHelpers::dimensions;
-using Opm::UgGridHelpers::cartDims;
-using Opm::UgGridHelpers::globalCell;
-using Opm::UgGridHelpers::cell2Faces;
-using Opm::UgGridHelpers::increment;
-using Opm::UgGridHelpers::getCoordinate;
+} //end namespace UgGridHelpers
+namespace AutoDiffGrid
+{
 
 /// \brief Mapps a grid type to the corresponding face to cell mapping.
 ///
@@ -54,16 +50,6 @@ struct ADFaceCellTraits
 {
 };
 
-template<>
-struct ADFaceCellTraits<UnstructuredGrid>
-{
-    typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> Type;
-};
-
-/// \brief Get the face to cell mapping of a grid.
-ADFaceCellTraits<UnstructuredGrid>::Type
-faceCells(const UnstructuredGrid& grid);
-
 /// \brief Get the z coordinates of the cell centroids of a grid.
 Eigen::Array<double, Eigen::Dynamic, 1>
 cellCentroidsZ(const UnstructuredGrid& grid);
@@ -120,10 +106,8 @@ void extractInternalFaces(const UnstructuredGrid& grid,
                           Eigen::Array<int, Eigen::Dynamic, 1>& internal_faces,
                           Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>& nbi);
 
-using Opm::UgGridHelpers::beginFaceCentroids;
-using Opm::UgGridHelpers::beginCellCentroids;
-}
-}
+} // end namespace AutoDiffGrid
+} // end namespace Opm
 
 #ifdef HAVE_DUNE_CORNERPOINT
 
@@ -133,24 +117,6 @@ namespace Opm
 
 namespace AutoDiffGrid
 {
-/// \brief Get the number of cells of a grid.
-int numCells(const Dune::CpGrid& grid);
-
-/// \brief Get the number of faces of a grid.
-int numFaces(const  Dune::CpGrid& grid);
-
-/// \brief Get the dimensions of a grid
-int dimensions(const Dune::CpGrid& grid);
-
-/// \brief Get the cartesion dimension of the underlying structured grid.
-const int* cartDims(const Dune::CpGrid& grid);
-
-/// \brief Get the local to global index mapping.
-///
-/// The global index is the index of the active cell
-/// in the underlying structured grid.
-const int*  globalCell(const Dune::CpGrid&);
-
 /// \brief A proxy class representing a row of FaceCellsContainer.
 class FaceCellsProxy
 {    
@@ -176,6 +142,8 @@ private:
 class FaceCellsContainerProxy
 {
 public:
+    typedef FaceCellsProxy row_type;
+    
     /// \brief Constructor.
     /// \param grid The grid whose information we represent.
     FaceCellsContainerProxy(const Dune::CpGrid* grid)
@@ -200,34 +168,11 @@ private:
     const Dune::CpGrid* grid_;
 };
 
-template<>
-struct ADFaceCellTraits<Dune::CpGrid>
-{
-    typedef FaceCellsContainerProxy Type;
-};
-
-/// \brief Get the face to cell mapping of a grid.
-ADFaceCellTraits<Dune::CpGrid>::Type
-faceCells(const Dune::CpGrid& grid);
-
-/// \brief Get the z coordinates of the cell centroids of a grid.
-Eigen::Array<double, Eigen::Dynamic, 1>
-cellCentroidsZ(const  Dune::CpGrid& grid);
-
-/// \brief Get the centroid of a cell.
-/// \param grid The grid whose cell centroid we query.
-/// \param cell_index The index of the corresponding cell.
-const double* cellCentroid(const Dune::CpGrid& grid, int cell_index);
-
-/// \brief Get the cell centroid of a face.
-/// \param grid The grid whose cell centroid we query.
-/// \param face_index The index of the corresponding face.
-const double* faceCentroid(const Dune::CpGrid& grid, int face_index);
-
 class Cell2FacesRow
 {
 public:
     class iterator
+        : public Dune::RandomAccessIteratorFacade<iterator,int, int, int>
     {
     public:
         iterator(const Dune::cpgrid::OrientedEntityTable<0,1>::row_type* row,
@@ -235,21 +180,35 @@ public:
             : row_(row), index_(index)
         {}
 
-        iterator operator++()
+        void increment()
         {
             ++index_;
-            return *this;
         }
-        iterator operator++(int)
+        void decrement()
         {
-            iterator ret=*this;
-            ++index_;
-            return ret;
+            --index_;
         }
-        int operator*()
+        int dereference() const
         {
             return row_->operator[](index_).index();
         }
+        int elementAt(int n) const
+        {
+            return row_->operator[](n).index();
+        }
+        void advance(int n)
+        {
+            index_+=n;
+        }
+        int distanceTo(const iterator& o)const
+        {
+            return o.index_-index_;
+        }
+        bool equals(const iterator& o) const
+        {
+            return index_==o.index_;
+        }
+        
     private:
         const Dune::cpgrid::OrientedEntityTable<0,1>::row_type* row_;
         int index_;
@@ -278,18 +237,169 @@ private:
 class Cell2FacesContainer
 {
 public:
+    typedef  Cell2FacesRow row_type;
+    
     Cell2FacesContainer(const Dune::CpGrid* grid)
         : grid_(grid)
     {};
     
-    Cell2FacesRow operator[](int cell_index)
+    Cell2FacesRow operator[](int cell_index) const
     {
         return Cell2FacesRow(grid_->cellFaceRow(cell_index));
     }
     
+        /// \brief Get the number of non-zero entries.
+    std::size_t noEntries() const
+    {
+        return grid_->numCellFaces();
+    }
 private:
     const Dune::CpGrid* grid_;
 };
+}
+
+namespace UgGridHelpers
+{
+template<>
+struct Cell2FacesTraits<Dune::CpGrid>
+{
+    typedef Opm::AutoDiffGrid::Cell2FacesContainer Type;
+};
+/// \brief An iterator over the cell volumes.
+template<const Dune::FieldVector<double, 3>& (Dune::CpGrid::*Method)(int)const>
+class CpGridCentroidIterator
+    : public Dune::RandomAccessIteratorFacade<CpGridCentroidIterator<Method>, Dune::FieldVector<double, 3>,
+                                              const Dune::FieldVector<double, 3>&, int>
+{
+public:
+    /// \brief Creates an iterator.
+    /// \param grid The grid the iterator belongs to.
+    /// \param cell_index The position of the iterator.
+    CpGridCentroidIterator(const  Dune::CpGrid& grid, int cell_index)
+        : grid_(&grid), cell_index_(cell_index)
+    {}
+
+    const Dune::FieldVector<double, 3>& dereference() const
+    {
+        return std::mem_fn(Method)(*grid_, cell_index_);
+    }
+    void increment()
+    {
+        ++cell_index_;
+    }
+    const Dune::FieldVector<double, 3>& elementAt(int n) const
+    {
+        return  std::mem_fn(Method)(*grid_, cell_index_);
+    }
+    void advance(int n)
+    {
+        cell_index_+=n;
+    }
+    void decrement()
+    {
+        --cell_index_;
+    }
+    int distanceTo(const CpGridCentroidIterator& o) const
+    {
+        return o.cell_index_-cell_index_;
+    }
+    bool equals(const CpGridCentroidIterator& o) const
+    {
+        return o.grid_==grid_ && o.cell_index_==cell_index_;
+    }
+    
+private:
+    const Dune::CpGrid* grid_;
+    int cell_index_;
+};
+
+template<>
+struct CellCentroidTraits<Dune::CpGrid>
+{
+    typedef CpGridCentroidIterator<&Dune::CpGrid::cellCentroid> IteratorType;
+    typedef const double* ValueType;
+};
+
+/// \brief Get the number of cells of a grid.
+int numCells(const Dune::CpGrid& grid);
+
+/// \brief Get the number of faces of a grid.
+int numFaces(const  Dune::CpGrid& grid);
+
+/// \brief Get the dimensions of a grid
+int dimensions(const Dune::CpGrid& grid);
+
+/// \brief Get the number of faces, where each face counts as many times as there are adjacent faces
+int numCellFaces(const Dune::CpGrid& grid);
+
+/// \brief Get the cartesion dimension of the underlying structured grid.
+const int* cartDims(const Dune::CpGrid& grid);
+
+/// \brief Get the local to global index mapping.
+///
+/// The global index is the index of the active cell
+/// in the underlying structured grid.
+const int*  globalCell(const Dune::CpGrid&);
+
+CellCentroidTraits<Dune::CpGrid>::IteratorType
+beginCellCentroids(const Dune::CpGrid& grid);
+
+/// \brief Get a coordinate of a specific cell centroid.
+/// \brief grid The grid.
+/// \brief cell_index The index of the specific cell.
+/// \breif coordinate The coordinate index.
+double cellCentroidCoordinate(const UnstructuredGrid& grid, int cell_index,
+                                 int coordinate);
+
+template<>
+struct FaceCentroidTraits<Dune::CpGrid>
+{
+    typedef CpGridCentroidIterator<&Dune::CpGrid::faceCentroid> IteratorType;
+    typedef const Dune::CpGrid::Vector ValueType;
+};
+
+/// \brief Get an iterator over the face centroids positioned at the first cell.
+FaceCentroidTraits<Dune::CpGrid>::IteratorType
+beginFaceCentroids(const Dune::CpGrid& grid);
+
+/// \brief Get a coordinate of a specific face centroid.
+/// \param grid The grid.
+/// \param face_index The index of the specific face.
+/// \param coordinate The coordinate index.
+FaceCentroidTraits<Dune::CpGrid>::ValueType
+faceCentroid(const Dune::CpGrid& grid, int face_index);
+
+template<>
+struct FaceCellTraits<Dune::CpGrid>
+{
+    typedef Opm::AutoDiffGrid::FaceCellsContainerProxy Type;
+};
+/// \brief Get the cell to faces mapping of a grid.
+Opm::AutoDiffGrid::Cell2FacesContainer cell2Faces(const Dune::CpGrid& grid);
+
+/// \brief Get the face to cell mapping of a grid.
+FaceCellTraits<Dune::CpGrid>::Type
+faceCells(const Dune::CpGrid& grid);
+
+const double* faceNormal(const Dune::CpGrid& grid, int face_index);
+} // end namespace UgGridHelperHelpers
+
+namespace AutoDiffGrid
+{
+
+/// \brief Get the z coordinates of the cell centroids of a grid.
+Eigen::Array<double, Eigen::Dynamic, 1>
+cellCentroidsZ(const  Dune::CpGrid& grid);
+
+/// \brief Get the centroid of a cell.
+/// \param grid The grid whose cell centroid we query.
+/// \param cell_index The index of the corresponding cell.
+const double* cellCentroid(const Dune::CpGrid& grid, int cell_index);
+
+/// \brief Get the cell centroid of a face.
+/// \param grid The grid whose cell centroid we query.
+/// \param face_index The index of the corresponding face.
+const double* faceCentroid(const Dune::CpGrid& grid, int face_index);
 
 template<>
 struct ADCell2FacesTraits<Dune::CpGrid>
@@ -297,9 +407,6 @@ struct ADCell2FacesTraits<Dune::CpGrid>
     typedef Cell2FacesContainer Type;
 };
 
-/// \brief Get the cell to faces mapping of a grid.
-Cell2FacesContainer cell2Faces(const Dune::CpGrid& grid);
-
 /// \brief Get the volume of a cell.
 /// \param grid The grid the cell belongs to.
 /// \param cell_index The index of the cell.
@@ -317,7 +424,7 @@ public:
         : grid_(&grid), cell_index_(cell_index)
     {}
 
-    double dereference()
+    double dereference() const
     {
         return grid_->cellVolume(cell_index_);
     }
@@ -325,11 +432,11 @@ public:
     {
         ++cell_index_;
     }
-    double elementAt(int n)
+    double elementAt(int n) const
     {
         return grid_->cellVolume(n);
     }
-    void adavance(int n)
+    void advance(int n)
     {
         cell_index_+=n;
     }
@@ -337,11 +444,11 @@ public:
     {
         --cell_index_;
     }
-    int distanceTo(const CellVolumeIterator& o)
+    int distanceTo(const CellVolumeIterator& o) const
     {
         return o.cell_index_-cell_index_;
     }
-    bool equals(const CellVolumeIterator& o)
+    bool equals(const CellVolumeIterator& o) const
     {
         return o.grid_==grid_ && o.cell_index_==cell_index_;
     }
@@ -362,8 +469,58 @@ CellVolumeIterator beginCellVolumes(const Dune::CpGrid& grid);
 
 /// \brief Get an iterator over the cell volumes of a grid positioned one after the last cell.
 CellVolumeIterator endCellVolumes(const Dune::CpGrid& grid);
+
+/// \brief extracts the internal faces of a grid.
+/// \param[in] The grid whose internal faces we query.
+/// \param[out] internal_faces The internal faces.
+/// \param[out] nbi 
+void extractInternalFaces(const Dune::CpGrid& grid,
+                          Eigen::Array<int, Eigen::Dynamic, 1>& internal_faces,
+                          Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>& nbi);
+
+template<>
+struct ADFaceCellTraits<Dune::CpGrid>
+    : public Opm::UgGridHelpers::FaceCellTraits<Dune::CpGrid>
+{};
+/// \brief Get the face to cell mapping of a grid.
+inline ADFaceCellTraits<Dune::CpGrid>::Type
+faceCells(const Dune::CpGrid& grid)
+{
+    return Opm::UgGridHelpers::faceCells(grid);
+}
 } // end namespace AutoDiffGrid
 } //end namespace OPM
 
 #endif
+namespace Opm
+{
+namespace AutoDiffGrid
+{
+
+using Opm::UgGridHelpers::SparseTableView;
+using Opm::UgGridHelpers::numCells;
+using Opm::UgGridHelpers::numFaces;
+using Opm::UgGridHelpers::dimensions;
+using Opm::UgGridHelpers::cartDims;
+using Opm::UgGridHelpers::globalCell;
+using Opm::UgGridHelpers::cell2Faces;
+using Opm::UgGridHelpers::increment;
+using Opm::UgGridHelpers::getCoordinate;
+using Opm::UgGridHelpers::numCellFaces;
+using Opm::UgGridHelpers::beginFaceCentroids;
+using Opm::UgGridHelpers::beginCellCentroids;
+
+template<>
+struct ADFaceCellTraits<UnstructuredGrid>
+{
+    typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> Type;
+};
+
+/// \brief Get the face to cell mapping of a grid.
+ADFaceCellTraits<UnstructuredGrid>::Type
+faceCells(const UnstructuredGrid& grid);
+
+}
+}
+
 #endif

opm/autodiff/ImpesTPFAAD.cpp

@@ -16,6 +16,7 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
+#include <config.h>
 
 #include <opm/autodiff/ImpesTPFAAD.hpp>
 #include <opm/autodiff/GeoProps.hpp>
@@ -61,11 +62,14 @@ namespace {
         std::vector<int> f2hf(2 * numFaces(grid), -1);
         Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
             face_cells;
-        SparseTableView c2f=cell2Faces(grid);
+        
+        typedef typename Opm::UgGridHelpers::Cell2FacesTraits<UnstructuredGrid>::Type
+            Cell2Faces;
+        Cell2Faces c2f=cell2Faces(grid);
         for (int c = 0; c < nc; ++c) {
-            typename SparseTableView::row_type 
+            typename Cell2Faces::row_type 
                 cell_faces = c2f[c];
-            typedef typename SparseTableView::row_type::iterator Iter;
+            typedef typename Cell2Faces::row_type::iterator Iter;
             for (Iter f=cell_faces.begin(), end=cell_faces.end();
                  f!=end; ++end) {
                 const int p = 0 + (face_cells(*f,0) != c);

opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -41,9 +41,12 @@ namespace Opm
     struct SimulatorReport;
 
     /// Class collecting all necessary components for a two-phase simulation.
+    template<class T>
     class SimulatorFullyImplicitBlackoil
     {
     public:
+        /// \brief The type of the grid that we use.
+        typedef T Grid;
         /// Initialise from parameters and objects to observe.
         /// \param[in] param       parameters, this class accepts the following:
         ///     parameter (default)            effect
@@ -67,7 +70,7 @@ namespace Opm
         /// \param[in] linsolver     linear solver
         /// \param[in] gravity       if non-null, gravity vector
         SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
-                                       const UnstructuredGrid& grid,
+                                       const Grid& grid,
                                        BlackoilPropsAdInterface& props,
                                        const RockCompressibility* rock_comp_props,
                                        WellsManager& wells_manager,
@@ -94,4 +97,5 @@ namespace Opm
 
 } // namespace Opm
 
+#include "SimulatorFullyImplicitBlackoil_impl.hpp"
 #endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp

@@ -0,0 +1,174 @@
+#include "config.h"
+
+#include "SimulatorFullyImplicitBlackoilOutput.hpp"
+
+#include <opm/core/utility/DataMap.hpp>
+#include <opm/core/io/vtk/writeVtkData.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+
+#include <opm/autodiff/GridHelpers.hpp>
+
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+
+#include <boost/filesystem.hpp>
+
+namespace Opm
+{
+
+
+    void outputStateVtk(const UnstructuredGrid& grid,
+                        const Opm::BlackoilState& 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(AutoDiffGrid::numCells(grid),
+                                  AutoDiffGrid::numFaces(grid),
+                                  AutoDiffGrid::beginFaceCentroids(grid),
+                                  AutoDiffGrid::faceCells(grid),
+                                  AutoDiffGrid::beginCellCentroids(grid),
+                                  AutoDiffGrid::beginCellVolumes(grid),
+                                  AutoDiffGrid::dimensions(grid),
+                                  state.faceflux(), cell_velocity);
+        dm["velocity"] = &cell_velocity;
+        Opm::writeVtkData(grid, dm, vtkfile);
+    }
+
+
+    void outputStateMatlab(const UnstructuredGrid& grid,
+                           const Opm::BlackoilState& state,
+                           const int step,
+                           const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["saturation"] = &state.saturation();
+        dm["pressure"] = &state.pressure();
+        dm["surfvolume"] = &state.surfacevol();
+        std::vector<double> cell_velocity;
+        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
+                                  AutoDiffGrid::numFaces(grid),
+                                  AutoDiffGrid::beginFaceCentroids(grid),
+                                  AutoDiffGrid::faceCells(grid),
+                                  AutoDiffGrid::beginCellCentroids(grid),
+                                  AutoDiffGrid::beginCellVolumes(grid),
+                                  AutoDiffGrid::dimensions(grid),
+                                  state.faceflux(), cell_velocity);
+        dm["velocity"] = &cell_velocity;
+
+        // Write data (not grid) in Matlab format
+        for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
+            std::ostringstream fname;
+            fname << output_dir << "/" << it->first;
+            boost::filesystem::path fpath = fname.str();
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+            }
+            fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+            std::ofstream file(fname.str().c_str());
+            if (!file) {
+                OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
+            }
+            file.precision(15);
+            const std::vector<double>& d = *(it->second);
+            std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
+        }
+    }
+    void outputWellStateMatlab(const Opm::WellState& well_state,
+                               const int step,
+                               const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["bhp"] = &well_state.bhp();
+        dm["wellrates"] = &well_state.wellRates();
+
+        // Write data (not grid) in Matlab format
+        for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
+            std::ostringstream fname;
+            fname << output_dir << "/" << it->first;
+            boost::filesystem::path fpath = fname.str();
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath);
+            }
+            fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+            std::ofstream file(fname.str().c_str());
+            if (!file) {
+                OPM_THROW(std::runtime_error,"Failed to open " << fname.str());
+            }
+            file.precision(15);
+            const std::vector<double>& d = *(it->second);
+            std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
+        }
+    }
+
+#if 0
+    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);
+    }
+
+    void outputWellReport(const Opm::WellReport& wellreport,
+                          const std::string& output_dir)
+    {
+        // Write well report.
+        std::string fname = output_dir  + "/wellreport.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        wellreport.write(os);
+    }
+#endif
+
+#ifdef HAVE_DUNE_CORNERPOINT
+    void outputStateVtk(const Dune::CpGrid& grid,
+                        const Opm::BlackoilState& state,
+                        const int step,
+                        const std::string& output_dir)
+    {
+        OPM_THROW(std::runtime_error, "outputStateVtk not implemented");
+    }
+    void outputStateMatlab(const Dune::CpGrid& grid,
+                           const Opm::BlackoilState& state,
+                           const int step,
+                           const std::string& output_dir)
+    {
+        OPM_THROW(std::runtime_error, "outputStateMatlab not implemented");
+    }
+
+#endif
+
+}

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -0,0 +1,42 @@
+#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED
+#include <opm/core/grid.h>
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+
+#include <string>
+
+#ifdef HAVE_DUNE_CORNERPOINT
+#include <dune/grid/CpGrid.hpp>
+#endif
+namespace Opm
+{
+
+    void outputStateVtk(const UnstructuredGrid& grid,
+                        const Opm::BlackoilState& state,
+                        const int step,
+                        const std::string& output_dir);
+
+
+    void outputStateMatlab(const UnstructuredGrid& grid,
+                           const Opm::BlackoilState& state,
+                           const int step,
+                           const std::string& output_dir);
+
+    void outputWellStateMatlab(const Opm::WellState& well_state,
+                               const int step,
+                               const std::string& output_dir);
+#ifdef HAVE_DUNE_CORNERPOINT
+    void outputStateVtk(const Dune::CpGrid& grid,
+                        const Opm::BlackoilState& state,
+                        const int step,
+                        const std::string& output_dir);
+
+
+    void outputStateMatlab(const Dune::CpGrid& grid,
+                           const Opm::BlackoilState& state,
+                           const int step,
+                           const std::string& output_dir);
+#endif
+}
+#endif

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -17,11 +17,6 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-
-#if HAVE_CONFIG_H
-#include "config.h"
-#endif // HAVE_CONFIG_H
-
 #include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
@@ -29,7 +24,7 @@
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
-
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
 #include <opm/core/pressure/flow_bc.h>
@@ -38,7 +33,6 @@
 #include <opm/core/simulator/SimulatorTimer.hpp>
 #include <opm/core/utility/StopWatch.hpp>
 #include <opm/core/io/eclipse/EclipseWriter.hpp>
-#include <opm/core/io/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 
@@ -62,12 +56,12 @@
 
 namespace Opm
 {
-
-    class SimulatorFullyImplicitBlackoil::Impl
+    template<class T>
+    class SimulatorFullyImplicitBlackoil<T>::Impl
     {
     public:
         Impl(const parameter::ParameterGroup& param,
-             const UnstructuredGrid& grid,
+             const Grid& grid,
              BlackoilPropsAdInterface& props,
              const RockCompressibility* rock_comp_props,
              WellsManager& wells_manager,
@@ -91,7 +85,7 @@ namespace Opm
         bool check_well_controls_;
         int max_well_control_iterations_;
         // Observed objects.
-        const UnstructuredGrid& grid_;
+        const Grid& grid_;
         BlackoilPropsAdInterface& props_;
         const RockCompressibility* rock_comp_props_;
         WellsManager& wells_manager_;
@@ -99,7 +93,7 @@ namespace Opm
         const double* gravity_;
         // Solvers
         DerivedGeology geo_;
-        FullyImplicitBlackoilSolver solver_;
+        FullyImplicitBlackoilSolver<Grid> solver_;
         // Misc. data
         std::vector<int> allcells_;
         EclipseWriter &eclipseWriter_;
@@ -108,8 +102,9 @@ namespace Opm
 
 
 
-    SimulatorFullyImplicitBlackoil::SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
-                                                                   const UnstructuredGrid& grid,
+    template<class T>
+    SimulatorFullyImplicitBlackoil<T>::SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
+                                                                   const Grid& grid,
                                                                    BlackoilPropsAdInterface& props,
                                                                    const RockCompressibility* rock_comp_props,
                                                                    WellsManager& wells_manager,
@@ -125,7 +120,8 @@ namespace Opm
 
 
 
-    SimulatorReport SimulatorFullyImplicitBlackoil::run(SimulatorTimer& timer,
+    template<class T>
+    SimulatorReport SimulatorFullyImplicitBlackoil<T>::run(SimulatorTimer& timer,
                                                         BlackoilState& state,
                                                         WellState& well_state)
     {
@@ -134,144 +130,11 @@ namespace Opm
 
 
 
-    static void outputStateVtk(const UnstructuredGrid& grid,
-                               const Opm::BlackoilState& 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(AutoDiffGrid::numCells(grid),
-                                  AutoDiffGrid::numFaces(grid),
-                                  AutoDiffGrid::beginFaceCentroids(grid),
-                                  AutoDiffGrid::faceCells(grid),
-                                  AutoDiffGrid::beginCellCentroids(grid),
-                                  AutoDiffGrid::beginCellVolumes(grid),
-                                  AutoDiffGrid::dimensions(grid),
-                                  state.faceflux(), cell_velocity);
-        dm["velocity"] = &cell_velocity;
-        Opm::writeVtkData(grid, dm, vtkfile);
-    }
-
-
-    static void outputStateMatlab(const UnstructuredGrid& grid,
-                                  const Opm::BlackoilState& state,
-                                  const int step,
-                                  const std::string& output_dir)
-    {
-        Opm::DataMap dm;
-        dm["saturation"] = &state.saturation();
-        dm["pressure"] = &state.pressure();
-        dm["surfvolume"] = &state.surfacevol();
-        std::vector<double> cell_velocity;
-        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
-                                  AutoDiffGrid::numFaces(grid),
-                                  AutoDiffGrid::beginFaceCentroids(grid),
-                                  AutoDiffGrid::faceCells(grid),
-                                  AutoDiffGrid::beginCellCentroids(grid),
-                                  AutoDiffGrid::beginCellVolumes(grid),
-                                  AutoDiffGrid::dimensions(grid),
-                                  state.faceflux(), cell_velocity);
-        dm["velocity"] = &cell_velocity;
-
-        // Write data (not grid) in Matlab format
-        for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
-            std::ostringstream fname;
-            fname << output_dir << "/" << it->first;
-            boost::filesystem::path fpath = fname.str();
-            try {
-                create_directories(fpath);
-            }
-            catch (...) {
-                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
-            }
-            fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
-            std::ofstream file(fname.str().c_str());
-            if (!file) {
-                OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
-            }
-            file.precision(15);
-            const std::vector<double>& d = *(it->second);
-            std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
-        }
-    }
-    static void outputWellStateMatlab(const Opm::WellState& well_state,
-                                  const int step,
-                                  const std::string& output_dir)
-    {
-        Opm::DataMap dm;
-        dm["bhp"] = &well_state.bhp();
-        dm["wellrates"] = &well_state.wellRates();
-
-        // Write data (not grid) in Matlab format
-        for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
-            std::ostringstream fname;
-            fname << output_dir << "/" << it->first;
-            boost::filesystem::path fpath = fname.str();
-            try {
-                create_directories(fpath);
-            }
-            catch (...) {
-                OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath);
-            }
-            fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
-            std::ofstream file(fname.str().c_str());
-            if (!file) {
-                OPM_THROW(std::runtime_error,"Failed to open " << fname.str());
-            }
-            file.precision(15);
-            const std::vector<double>& d = *(it->second);
-            std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
-        }
-    }
-
-#if 0
-    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);
-    }
-
-    static void outputWellReport(const Opm::WellReport& wellreport,
-                                 const std::string& output_dir)
-    {
-        // Write well report.
-        std::string fname = output_dir  + "/wellreport.txt";
-        std::ofstream os(fname.c_str());
-        if (!os) {
-            OPM_THROW(std::runtime_error, "Failed to open " << fname);
-        }
-        wellreport.write(os);
-    }
-#endif
-
 
     // \TODO: Treat bcs.
-    SimulatorFullyImplicitBlackoil::Impl::Impl(const parameter::ParameterGroup& param,
-                                               const UnstructuredGrid& grid,
+    template<class T>
+    SimulatorFullyImplicitBlackoil<T>::Impl::Impl(const parameter::ParameterGroup& param,
+                                               const Grid& grid,
                                                BlackoilPropsAdInterface& props,
                                                const RockCompressibility* rock_comp_props,
                                                WellsManager& wells_manager,
@@ -321,10 +184,8 @@ namespace Opm
         }
     }
 
-
-
-
-    SimulatorReport SimulatorFullyImplicitBlackoil::Impl::run(SimulatorTimer& timer,
+    template<class T>
+    SimulatorReport SimulatorFullyImplicitBlackoil<T>::Impl::run(SimulatorTimer& timer,
                                                               BlackoilState& state,
                                                               WellState& well_state)
     {
@@ -334,9 +195,9 @@ namespace Opm
         // Initialisation.
         std::vector<double> porevol;
         if (rock_comp_props_ && rock_comp_props_->isActive()) {
-            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+            computePorevolume(AutoDiffGrid::numCells(grid_), AutoDiffGrid::beginCellVolumes(grid_), props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
         } else {
-            computePorevolume(grid_, props_.porosity(), porevol);
+            computePorevolume(AutoDiffGrid::numCells(grid_), AutoDiffGrid::beginCellVolumes(grid_), props_.porosity(), porevol);
         }
         // const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
         std::vector<double> initial_porevol = porevol;
@@ -420,7 +281,7 @@ namespace Opm
             // Update pore volumes if rock is compressible.
             if (rock_comp_props_ && rock_comp_props_->isActive()) {
                 initial_porevol = porevol;
-                computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+                computePorevolume(AutoDiffGrid::numCells(grid_), AutoDiffGrid::beginCellVolumes(grid_), props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
             }
 
             // Hysteresis