Merge pull request #100 from blattms/master-refactor-for-cpgrid-support

Adds fully implicite black oil solver working with CpGrid

CMakeLists.txt

@@ -46,13 +46,18 @@ include (${project}-prereqs)
 include (CMakeLists_files.cmake)
 
 macro (config_hook)
-#	opm_need_version_of ("dune-common")
+	opm_need_version_of ("dune-common")
 endmacro (config_hook)
 
 macro (prereqs_hook)
 endmacro (prereqs_hook)
 
 macro (sources_hook)
+	if(DUNE_CORNERPOINT_FOUND OR dune-cornerpoint_FOUND)
+		list (APPEND examples_SOURCES
+			${PROJECT_SOURCE_DIR}/examples/sim_fibo_ad_cp.cpp
+			)
+	endif()
 endmacro (sources_hook)
 
 macro (fortran_hook)

CMakeLists_files.cmake

@@ -28,10 +28,10 @@
 list (APPEND MAIN_SOURCE_FILES
 	opm/autodiff/BlackoilPropsAd.cpp
 	opm/autodiff/BlackoilPropsAdInterface.cpp
-	opm/autodiff/FullyImplicitBlackoilSolver.cpp
+	opm/autodiff/GridHelpers.cpp
 	opm/autodiff/ImpesTPFAAD.cpp
 	opm/autodiff/SimulatorCompressibleAd.cpp
-	opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
+	opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
 	opm/autodiff/SimulatorIncompTwophaseAd.cpp
 	opm/autodiff/TransportSolverTwophaseAd.cpp
 	opm/autodiff/BlackoilPropsAdFromDeck.cpp
@@ -99,10 +99,13 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/BlackoilPropsAdFromDeck.hpp
 	opm/autodiff/BlackoilPropsAdInterface.hpp
 	opm/autodiff/GeoProps.hpp
+	opm/autodiff/GridHelpers.hpp
 	opm/autodiff/ImpesTPFAAD.hpp
 	opm/autodiff/FullyImplicitBlackoilSolver.hpp
+	opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
 	opm/autodiff/SimulatorCompressibleAd.hpp
 	opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
+	opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
 	opm/autodiff/SimulatorIncompTwophaseAd.hpp
 	opm/autodiff/TransportSolverTwophaseAd.hpp
 	opm/autodiff/WellDensitySegmented.hpp

cmake/Modules/Finddune-cornerpoint.cmake

@@ -26,7 +26,7 @@ find_opm_package (
   "dunecornerpoint"
 
   # defines to be added to compilations
-  ""
+  "HAVE_DUNE_CORNERPOINT"
 
   # test program
 "#include <dune/grid/CpGrid.hpp>

cmake/Modules/opm-autodiff-prereqs.cmake

@@ -3,6 +3,7 @@
 
 # defines that must be present in config.h for our headers
 set (opm-autodiff_CONFIG_VAR
+  HAVE_DUNE_CORNERPOINT
 	)
 
 # dependencies
@@ -17,6 +18,7 @@ set (opm-autodiff_DEPS
 	# DUNE prerequisites
 	"dune-common REQUIRED;
 	dune-istl REQUIRED;
+        dune-cornerpoint;
 	opm-core REQUIRED"
 	# Eigen
 	"Eigen3 3.1.2 "

dune.module

@@ -4,3 +4,4 @@ Version: 0.9
 Label: 2013.10
 Maintainer: atgeirr@sintef.no
 Depends: opm-core
+Suggests: dune-cornerpoint

examples/sim_fibo_ad.cpp

@@ -205,7 +205,7 @@ try
         }
 
         // Create and run simulator.
-        SimulatorFullyImplicitBlackoil simulator(param,
+        SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
                                                  *grid->c_grid(),
                                                  *new_props,
                                                  rock_comp->isActive() ? rock_comp.get() : 0,

examples/sim_fibo_ad_cp.cpp

@@ -0,0 +1,290 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+
+  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/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <dune/common/version.hh>
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
+#include <dune/common/parallel/mpihelper.hh>
+#else
+#include <dune/common/mpihelper.hh>
+#endif
+
+#include <opm/core/pressure/FlowBCManager.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/cornerpoint_grid.h>
+
+#include <opm/core/grid/GridManager.hpp>
+
+#include <dune/grid/CpGrid.hpp>
+#include <dune/grid/common/GridAdapter.hpp>
+
+#include <opm/core/wells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/simulator/initState.hpp>
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/core/io/eclipse/EclipseWriter.hpp>
+#include <opm/core/props/BlackoilPropertiesBasic.hpp>
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+
+#include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
+#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
+#include <opm/autodiff/GridHelpers.hpp>
+#include <opm/core/utility/share_obj.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/algorithm/string.hpp>
+
+#include <memory>
+#include <algorithm>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+#include <numeric>
+
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "----------------------------------------------------------------" << std::endl;
+        }
+    }
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    Dune::MPIHelper& helper= Dune::MPIHelper::instance(argc, argv);
+    using namespace Opm;
+
+    std::cout << "\n================    Test program for fully implicit three-phase black-oil flow     ===============\n\n";
+    parameter::ParameterGroup param(argc, argv, false);
+    std::cout << "---------------    Reading parameters     ---------------" << std::endl;
+
+    // If we have a "deck_filename", grid and props will be read from that.
+    bool use_deck = param.has("deck_filename");
+    if (!use_deck) {
+        OPM_THROW(std::runtime_error, "This program must be run with an input deck. "
+                  "Specify the deck with deck_filename=deckname.data (for example).");
+    }
+    std::shared_ptr<Dune::CpGrid> grid;
+    std::shared_ptr<BlackoilPropertiesInterface> props;
+    std::shared_ptr<BlackoilPropsAdInterface> new_props;
+    std::shared_ptr<RockCompressibility> rock_comp;
+    BlackoilState state;
+    // bool check_well_controls = false;
+    // int max_well_control_iterations = 0;
+    double gravity[3] = { 0.0 };
+    std::string deck_filename = param.get<std::string>("deck_filename");
+
+    Opm::ParserPtr newParser(new Opm::Parser() );
+    Opm::DeckConstPtr newParserDeck = newParser->parseFile( deck_filename );
+
+    // Grid init
+    grid.reset(new Dune::CpGrid());
+    {
+        grdecl g = {};
+        GridManager::createGrdecl(newParserDeck, g);
+
+        grid->processEclipseFormat(g, 2e-12, false);
+
+        std::free(const_cast<double*>(g.mapaxes));
+    }
+
+
+    Opm::EclipseWriter outputWriter(param, newParserDeck,
+                                    Opm::UgGridHelpers::numCells(*grid),
+                                    Opm::UgGridHelpers::globalCell(*grid),
+                                    Opm::UgGridHelpers::cartDims(*grid),
+                                    Opm::UgGridHelpers::dimensions(*grid));
+
+    // Rock and fluid init
+    props.reset(new BlackoilPropertiesFromDeck(newParserDeck, Opm::UgGridHelpers::numCells(*grid),
+                                               Opm::UgGridHelpers::globalCell(*grid),
+                                               Opm::UgGridHelpers::cartDims(*grid),
+                                               Opm::UgGridHelpers::beginCellCentroids(*grid),
+                                               Opm::UgGridHelpers::dimensions(*grid), param));
+    new_props.reset(new BlackoilPropsAdFromDeck(newParserDeck, *grid));
+    // check_well_controls = param.getDefault("check_well_controls", false);
+    // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+    // Rock compressibility.
+    rock_comp.reset(new RockCompressibility(newParserDeck));
+
+    // Gravity.
+    gravity[2] = newParserDeck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
+
+    // Init state variables (saturation and pressure).
+    if (param.has("init_saturation")) {
+        initStateBasic(grid->numCells(), &(grid->globalCell())[0],
+                       &(grid->logicalCartesianSize()[0]),
+                       grid->numFaces(), AutoDiffGrid::faceCells(*grid),
+                       grid->beginFaceCentroids(),
+                       grid->beginCellCentroids(), Dune::CpGrid::dimension,
+                       *props, param, gravity[2], state);
+        initBlackoilSurfvol(grid->numCells(), *props, state);
+        enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
+        const PhaseUsage pu = props->phaseUsage();
+        if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
+            const int np = props->numPhases();
+            const int nc = grid->numCells();
+            for (int c = 0; c < nc; ++c) {
+                state.gasoilratio()[c] = state.surfacevol()[c*np + pu.phase_pos[Gas]]
+                    / state.surfacevol()[c*np + pu.phase_pos[Oil]];
+            }
+        }
+    } else {
+        initBlackoilStateFromDeck(grid->numCells(), &(grid->globalCell())[0],
+                                  grid->numFaces(), AutoDiffGrid::faceCells(*grid),
+                                  grid->beginFaceCentroids(),
+                                  grid->beginCellCentroids(), Dune::CpGrid::dimension,
+                                  *props, newParserDeck, gravity[2], state);
+    }
+
+    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
+    const double *grav = use_gravity ? &gravity[0] : 0;
+
+    // Linear solver.
+    LinearSolverFactory linsolver(param);
+
+    // Write parameters used for later reference.
+    bool output = param.getDefault("output", true);
+    std::ofstream outStream;
+    std::string output_dir;
+    if (output) {
+        output_dir =
+            param.getDefault("output_dir", std::string("output"));
+        boost::filesystem::path fpath(output_dir);
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        std::string filename = output_dir + "/timing.param";
+        outStream.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
+        // open file to clean it. The file is appended to in SimulatorTwophase
+        filename = output_dir + "/step_timing.param";
+        std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
+        step_os.close();
+        param.writeParam(output_dir + "/simulation.param");
+    }
+
+    std::cout << "\n\n================    Starting main simulation loop     ===============\n"
+              << std::flush;
+
+    WellStateFullyImplicitBlackoil well_state;
+    Opm::TimeMapPtr timeMap(new Opm::TimeMap(newParserDeck));
+    SimulatorTimer simtimer;
+    std::shared_ptr<EclipseState> eclipseState(new EclipseState(newParserDeck));
+
+    // initialize variables
+    simtimer.init(timeMap);
+
+    SimulatorReport fullReport;
+    for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
+        // Report on start of a report step.
+        std::cout << "\n"
+                  << "---------------------------------------------------------------\n"
+                  << "--------------    Starting report step " << reportStepIdx << "    --------------\n"
+                  << "---------------------------------------------------------------\n"
+                  << "\n";
+
+        // Create new wells, well_state
+        WellsManager wells(eclipseState, reportStepIdx, Opm::UgGridHelpers::numCells(*grid), 
+                           Opm::UgGridHelpers::globalCell(*grid),
+                           Opm::UgGridHelpers::cartDims(*grid),
+                           Opm::UgGridHelpers::dimensions(*grid),
+                           Opm::UgGridHelpers::beginCellCentroids(*grid),
+                           Opm::UgGridHelpers::cell2Faces(*grid),
+                           Opm::UgGridHelpers::beginFaceCentroids(*grid),
+                           props->permeability());
+
+        if (reportStepIdx == 0) {
+            // @@@ HACK: we should really make a new well state and
+            // properly transfer old well state to it every epoch,
+            // since number of wells may change etc.
+            well_state.init(wells.c_wells(), state);
+        }
+
+        simtimer.setCurrentStepNum(reportStepIdx);
+
+        if (reportStepIdx == 0) {
+            outputWriter.writeInit(simtimer);
+            outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
+        }
+
+        // Create and run simulator.
+        SimulatorFullyImplicitBlackoil<Dune::CpGrid> simulator(param,
+                                                               *grid,
+                                                               *new_props,
+                                                               rock_comp->isActive() ? rock_comp.get() : 0,
+                                                               wells,
+                                                               linsolver,
+                                                               grav);
+        SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
+
+        ++simtimer;
+
+        outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
+        fullReport += episodeReport;
+    }
+
+    std::cout << "\n\n================    End of simulation     ===============\n\n";
+    fullReport.report(std::cout);
+
+    if (output) {
+        std::string filename = output_dir + "/walltime.param";
+        std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
+        fullReport.reportParam(tot_os);
+        warnIfUnusedParams(param);
+    }
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}
+

examples/test_implicit_ad.cpp

@@ -104,7 +104,7 @@ try
 
     Opm::LinearSolverFactory linsolver(param);
 
-    Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, 0, *wells, linsolver);
+    Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(*g, props, geo, 0, *wells, linsolver);
 
     Opm::BlackoilState state;
     initStateBasic(*g, props0, param, 0.0, state);

opm/autodiff/AutoDiffHelpers.hpp

@@ -21,6 +21,7 @@
 #define OPM_AUTODIFFHELPERS_HEADER_INCLUDED
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/GridHelpers.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/utility/ErrorMacros.hpp>
 
@@ -58,30 +59,20 @@ struct HelperOps
     M fulldiv;
 
     /// Constructs all helper vectors and matrices.
-    HelperOps(const UnstructuredGrid& grid)
+    template<class Grid>
+    HelperOps(const Grid& grid)
     {
-        const int nc = grid.number_of_cells;
+        using namespace AutoDiffGrid;
-        const int nf = grid.number_of_faces;
+        const int nc = numCells(grid);
+        const int nf = numFaces(grid);
         // 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);
+        TwoColInt nbi;
-        // std::cout << "nb = \n" << nb << std::endl;
+        extractInternalFaces(grid, internal_faces, nbi);
-        TwoColBool nbib = nb >= 0;
+        int num_internal=internal_faces.size();
-        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);
@@ -103,6 +94,7 @@ struct HelperOps
         div = ngrad.transpose();
         std::vector<Tri> fullngrad_tri;
         fullngrad_tri.reserve(2*nf);
+        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);
@@ -127,12 +119,16 @@ 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();
+            typename ADFaceCellTraits<Grid>::Type
+                face_cells = faceCells(g);
             assert(nif == ifaceflux.size());
 
             // Define selector structure.
@@ -140,8 +136,8 @@ struct HelperOps
             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 c1 = face_cells(f,0);
-                const int c2 = g.face_cells[2*f + 1];
+                const int c2 = face_cells(f,1);
 
                 assert ((c1 >= 0) && (c2 >= 0));
 
@@ -152,7 +148,7 @@ struct HelperOps
             }
 
             // Assemble explicit selector operator.
-            select_.resize(nif, g.number_of_cells);
+            select_.resize(nif, numCells(g));
             select_.setFromTriplets(s.begin(), s.end());
         }
 

opm/autodiff/BlackoilPropsAdFromDeck.cpp

@@ -51,9 +51,68 @@ namespace Opm
     BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
                                                      const UnstructuredGrid& grid,
                                                      const bool init_rock)
+    {
+        init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, 
+             grid.cell_centroids, grid.dimensions, init_rock);
+    }
+
+    /// Constructor wrapping an opm-core black oil interface.
+    BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
+                                                     const UnstructuredGrid& grid,
+                                                     const bool init_rock)
+    {
+        init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims, 
+             grid.cell_centroids, grid.dimensions, init_rock);
+    }
+
+#ifdef HAVE_DUNE_CORNERPOINT
+        
+    /// Constructor wrapping an opm-core black oil interface.
+    BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
+                                                     const Dune::CpGrid& grid,
+                                                     const bool init_rock )
+    {
+        init(deck, grid.numCells(), static_cast<const int*>(&grid.globalCell()[0]),
+             static_cast<const int*>(&grid.logicalCartesianSize()[0]),
+             grid.beginCellCentroids(), Dune::CpGrid::dimension, init_rock);
+    }
+
+    /// Constructor wrapping an opm-core black oil interface.
+    BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
+                                                     const Dune::CpGrid& grid,
+                                                     const bool init_rock )
+    {
+        init(newParserDeck, grid.numCells(), static_cast<const int*>(&grid.globalCell()[0]),
+             static_cast<const int*>(&grid.logicalCartesianSize()[0]),
+             grid.beginCellCentroids(), Dune::CpGrid::dimension, init_rock);
+    }
+#endif
+
+    /// Constructor wrapping an opm-core black oil interface.
+    template<class T>
+    BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
+                                                     int number_of_cells,
+                                                     const int* global_cell,
+                                                     const int* cart_dims,
+                                                     T begin_cell_centroids,
+                                                     int dimensions,
+                                                     const bool init_rock)
+    {
+        init(newParserDeck, number_of_cells, global_cell, cart_dims, begin_cell_centroids,
+             dimensions, init_rock);
+    }
+
+    template<class T>
+    void BlackoilPropsAdFromDeck::init(const EclipseGridParser& deck,
+                                       int number_of_cells,
+                                       const int* global_cell,
+                                       const int* cart_dims,
+                                       T begin_cell_centroids,
+                                       int dimensions,
+                                       const bool init_rock)
     {
         if (init_rock){
-            rock_.init(deck, grid);
+            rock_.init(deck, number_of_cells, global_cell, cart_dims);
         }
         const int samples = 0;
         const int region_number = 0;
@@ -122,7 +181,7 @@ namespace Opm
         SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
             = new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
         satprops_.reset(ptr);
-        ptr->init(deck, grid, -1);
+        ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions, -1);
 
         if (phase_usage_.num_phases != satprops_->numPhases()) {
             OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "
@@ -131,14 +190,17 @@ namespace Opm
         }
     }
 
-
+    template<class T>
-    /// Constructor wrapping an opm-core black oil interface.
+    void BlackoilPropsAdFromDeck::init(Opm::DeckConstPtr newParserDeck,
-    BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
+                                       int number_of_cells,
-                                                     const UnstructuredGrid& grid,
+                                       const int* global_cell,
-                                                     const bool init_rock)
+                                       const int* cart_dims,
+                                       T begin_cell_centroids,
+                                       int dimensions,
+                                       const bool init_rock)
     {
         if (init_rock){
-            rock_.init(newParserDeck, grid);
+            rock_.init(newParserDeck, number_of_cells, global_cell, cart_dims);
         }
 
         const int samples = 0;
@@ -218,7 +280,7 @@ namespace Opm
         SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
             = new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
         satprops_.reset(ptr);
-        ptr->init(newParserDeck, grid, -1);
+        ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids, dimensions, -1);
 
         if (phase_usage_.num_phases != satprops_->numPhases()) {
             OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "

opm/autodiff/BlackoilPropsAdFromDeck.hpp

@@ -33,6 +33,10 @@
 #include <boost/shared_ptr.hpp>
 #include <boost/scoped_ptr.hpp>
 
+#ifdef HAVE_DUNE_CORNERPOINT
+#include <dune/grid/CpGrid.hpp>
+#endif
+
 namespace Opm
 {
 
@@ -60,6 +64,41 @@ namespace Opm
                                 const UnstructuredGrid& grid,
                                 const bool init_rock = true );
 
+#ifdef HAVE_DUNE_CORNERPOINT
+        
+        /// Constructor wrapping an opm-core black oil interface.
+        BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
+                                const Dune::CpGrid& grid,
+                                const bool init_rock = true );
+
+        /// Constructor wrapping an opm-core black oil interface.
+        BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
+                                const Dune::CpGrid& grid,
+                                const bool init_rock = true );
+
+
+#endif
+
+        /// Constructor taking not a grid but only the needed information
+        template<class T>
+        BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
+                                int number_of_cells,
+                                const int* global_cell,
+                                const int* cart_dims,
+                                T begin_cell_centroids,
+                                int dimensions,
+                                const bool init_rock);
+
+        /// Constructor taking not a grid but only the needed information
+        template<class T>
+        BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
+                                int number_of_cells,
+                                const int* global_cell,
+                                const int* cart_dims,
+                                T begin_cell_centroids,
+                                int dimensions,
+                                const bool init_rock);
+
         ////////////////////////////
         //      Rock interface    //
         ////////////////////////////
@@ -329,6 +368,24 @@ namespace Opm
                            const std::vector<int>& cells);
 
     private:
+        /// Initializes the properties.
+        template<class T>
+        void init(const EclipseGridParser& deck,
+                  int number_of_cells,
+                  const int* global_cell,
+                  const int* cart_dims,
+                  T begin_cell_centroids,
+                  int dimension,
+                  const bool init_rock);
+        /// Initializes the properties.
+        template<class T>
+        void init(Opm::DeckConstPtr deck,
+                  int number_of_cells,
+                  const int* global_cell,
+                  const int* cart_dims,
+                  T begin_cell_centroids,
+                  int dimension,
+                  const bool init_rock);
         RockFromDeck rock_;
         boost::scoped_ptr<SaturationPropsInterface> satprops_;
         PhaseUsage phase_usage_;

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_;
@@ -269,5 +272,6 @@ namespace Opm {
     };
 } // namespace Opm
 
+#include "FullyImplicitBlackoilSolver_impl.hpp"
 
 #endif // OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -21,6 +21,7 @@
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/GridHelpers.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/WellDensitySegmented.hpp>
@@ -74,21 +75,27 @@ 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 )
     {
-        const int nc = grid.number_of_cells;
+        using namespace Opm::AutoDiffGrid;
+        const int nc = numCells(grid);
+        typedef typename Opm::UgGridHelpers::Cell2FacesTraits<Grid>::Type Cell2Faces;
+        Cell2Faces c2f = cell2Faces(grid);
 
-        std::vector<int> f2hf(2 * grid.number_of_faces, -1);
+        std::vector<int> f2hf(2 * numFaces(grid), -1);
+        typename ADFaceCellTraits<Grid>::Type
+            face_cells = faceCells(grid);
         for (int c = 0, i = 0; c < nc; ++c) {
-            for (; i < grid.cell_facepos[c + 1]; ++i) {
+            typename Cell2Faces::row_type faces=c2f[c];
-                const int f = grid.cell_faces[ i ];
+            typedef typename Cell2Faces::row_type::iterator Iter;
-                const int p = 0 + (grid.face_cells[2*f + 0] != c);
+            for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f) {
+                const int p = 0 + (face_cells(*f, 0) != c);
 
-                f2hf[2*f + p] = i;
+                f2hf[2*(*f) + p] = i;
             }
         }
 
@@ -104,8 +111,8 @@ namespace {
         std::vector<Tri> grav;  grav.reserve(2 * ni);
         for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
             const int f  = ops.internal_faces[ i ];
-            const int c1 = grid.face_cells[2*f + 0];
+            const int c1 = faceCells(grid)(f, 0);
-            const int c2 = grid.face_cells[2*f + 1];
+            const int c2 = faceCells(grid)(f, 1);
 
             assert ((c1 >= 0) && (c2 >= 0));
 
@@ -123,12 +130,13 @@ namespace {
     }
 
 
-
+    template<class Grid>
-    V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav)
+    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;
         const int nperf = wells.well_connpos[nw];
-        const int dim = grid.dimensions;
+        const int dim = dimensions(grid);
         V wdp = V::Zero(nperf,1);
         assert(wdp.size() == rho.size());
 
@@ -142,7 +150,7 @@ namespace {
             const double ref_depth = wells.depth_ref[w];
             for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) {
                 const int cell = wells.well_cells[j];
-                const double cell_depth = grid.cell_centroids[dim * cell + dim - 1];
+                const double cell_depth = cellCentroid(grid, cell)[dim - 1];
                 wdp[j] = rho[j]*grav*(cell_depth - ref_depth);
             }
         }
@@ -188,9 +196,9 @@ namespace {
 
 
 
-
+    template<class T>
-    FullyImplicitBlackoilSolver::
+    FullyImplicitBlackoilSolver<T>::
-    FullyImplicitBlackoilSolver(const UnstructuredGrid&         grid ,
+    FullyImplicitBlackoilSolver(const Grid&                     grid ,
                                 const BlackoilPropsAdInterface& fluid,
                                 const DerivedGeology&           geo  ,
                                 const RockCompressibility*      rock_comp_props,
@@ -204,12 +212,12 @@ namespace {
         , linsolver_ (linsolver)
         , active_(activePhases(fluid.phaseUsage()))
         , canph_ (active2Canonical(fluid.phaseUsage()))
-        , cells_ (buildAllCells(grid.number_of_cells))
+        , cells_ (buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
         , ops_   (grid)
         , 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() } )
@@ -219,9 +227,9 @@ namespace {
 
 
 
-
+    template<class T>
     void
-    FullyImplicitBlackoilSolver::
+    FullyImplicitBlackoilSolver<T>::
     step(const double   dt,
          BlackoilState& x ,
          WellStateFullyImplicitBlackoil& xw)
@@ -273,7 +281,8 @@ namespace {
 
 
 
-    FullyImplicitBlackoilSolver::ReservoirResidualQuant::ReservoirResidualQuant()
+    template<class T>
+    FullyImplicitBlackoilSolver<T>::ReservoirResidualQuant::ReservoirResidualQuant()
         : accum(2, ADB::null())
         , mflux(   ADB::null())
         , b    (   ADB::null())
@@ -286,7 +295,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())
@@ -300,7 +310,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)
@@ -331,9 +342,10 @@ namespace {
 
 
 
-    FullyImplicitBlackoilSolver::SolutionState
+    template<class T>
-    FullyImplicitBlackoilSolver::constantState(const BlackoilState& x,
+    typename FullyImplicitBlackoilSolver<T>::SolutionState
-                                               const WellStateFullyImplicitBlackoil& xw)
+    FullyImplicitBlackoilSolver<T>::constantState(const BlackoilState& x,
+                                                  const WellStateFullyImplicitBlackoil&     xw)
     {
         auto state = variableState(x, xw);
 
@@ -356,11 +368,13 @@ namespace {
 
 
 
-    FullyImplicitBlackoilSolver::SolutionState
+    template<class T>
-    FullyImplicitBlackoilSolver::variableState(const BlackoilState& x,
+    typename FullyImplicitBlackoilSolver<T>::SolutionState
-                                               const WellStateFullyImplicitBlackoil& xw)
+    FullyImplicitBlackoilSolver<T>::variableState(const BlackoilState& x,
+                                                  const WellStateFullyImplicitBlackoil&     xw)
     {
-        const int nc = grid_.number_of_cells;
+        using namespace Opm::AutoDiffGrid;
+        const int nc = numCells(grid_);
         const int np = x.numPhases();
 
         std::vector<V> vars0;
@@ -502,8 +516,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();
@@ -543,9 +558,11 @@ namespace {
 
 
 
-    void FullyImplicitBlackoilSolver::computeWellConnectionPressures(const SolutionState& state,
+    template<class T>
-                                                                     const WellStateFullyImplicitBlackoil& xw)
+    void FullyImplicitBlackoilSolver<T>::computeWellConnectionPressures(const SolutionState& state,
+                                                                        const WellStateFullyImplicitBlackoil& xw)
     {
+        using namespace Opm::AutoDiffGrid;
         // 1. Compute properties required by computeConnectionPressureDelta().
         //    Note that some of the complexity of this part is due to the function
         //    taking std::vector<double> arguments, and not Eigen objects.
@@ -583,7 +600,7 @@ namespace {
         // b is row major, so can just copy data.
         std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases);
         // Extract well connection depths.
-        const V depth = Eigen::Map<DataBlock>(grid_.cell_centroids, grid_.number_of_cells, grid_.dimensions).rightCols<1>();
+        const V depth = cellCentroidsZ(grid_);
         const V pdepth = subset(depth, well_cells);
         std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
         // Surface density.
@@ -591,7 +608,7 @@ namespace {
         // Gravity
         double grav = 0.0;
         const double* g = geo_.gravity();
-        const int dim = grid_.dimensions;
+        const int dim = dimensions(grid_);
         if (g) {
             // Guard against gravity in anything but last dimension.
             for (int dd = 0; dd < dim - 1; ++dd) {
@@ -612,12 +629,14 @@ namespace {
 
 
 
+    template<class T>
     void
-    FullyImplicitBlackoilSolver::
+    FullyImplicitBlackoilSolver<T>::
     assemble(const V&             pvdt,
              const BlackoilState& x   ,
              WellStateFullyImplicitBlackoil& xw  )
     {
+        using namespace Opm::AutoDiffGrid;
         // Create the primary variables.
         const SolutionState state = variableState(x, xw);
 
@@ -684,10 +703,11 @@ namespace {
 
 
 
-    void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state,
+    template <class T>
+    void FullyImplicitBlackoilSolver<T>::addWellEq(const SolutionState& state,
                                                 WellStateFullyImplicitBlackoil& xw)
     {
-        const int nc = grid_.number_of_cells;
+        const int nc = Opm::AutoDiffGrid::numCells(grid_);
         const int np = wells_.number_of_phases;
         const int nw = wells_.number_of_wells;
         const int nperf = wells_.well_connpos[nw];
@@ -935,8 +955,8 @@ namespace {
 
 
 
-
+    template<class T>
-    void FullyImplicitBlackoilSolver::updateWellControls(const ADB& bhp,
+    void FullyImplicitBlackoilSolver<T>::updateWellControls(const ADB& bhp,
                                                          const ADB& well_phase_flow_rate,
                                                          WellStateFullyImplicitBlackoil& xw) const
     {
@@ -989,8 +1009,8 @@ namespace {
 
 
 
-
+    template<class T>
-    void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state,
+    void FullyImplicitBlackoilSolver<T>::addWellControlEq(const SolutionState& state,
                                                        const WellStateFullyImplicitBlackoil& xw)
     {
         // Handling BHP and SURFACE_RATE wells.
@@ -1035,13 +1055,14 @@ namespace {
 
 
 
-    void FullyImplicitBlackoilSolver::addOldWellEq(const SolutionState& state)
+    template<class T>
+    void FullyImplicitBlackoilSolver<T>::addOldWellEq(const SolutionState& state)
     {
         // -------- Well equation, and well contributions to the mass balance equations --------
 
         // Contribution to mass balance will have to wait.
 
-        const int nc = grid_.number_of_cells;
+        const int nc = numCells(grid_);
         const int np = wells_.number_of_phases;
         const int nw = wells_.number_of_wells;
         const int nperf = wells_.well_connpos[nw];
@@ -1061,7 +1082,7 @@ namespace {
         // Compute well pressure differentials.
         // Construct pressure difference vector for wells.
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
-        const int dim = grid_.dimensions;
+        const int dim = dimensions(grid_);
         const double* g = geo_.gravity();
         if (g) {
             // Guard against gravity in anything but last dimension.
@@ -1156,7 +1177,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];
@@ -1208,12 +1230,14 @@ namespace {
 
 
 
-    void FullyImplicitBlackoilSolver::updateState(const V& dx,
+    template<class T>
+    void FullyImplicitBlackoilSolver<T>::updateState(const V& dx,
                                                   BlackoilState& state,
                                                   WellStateFullyImplicitBlackoil& well_state)
     {
+        using namespace Opm::AutoDiffGrid;
         const int np = fluid_.numPhases();
-        const int nc = grid_.number_of_cells;
+        const int nc = numCells(grid_);
         const int nw = wells_.number_of_wells;
         const V null;
         assert(null.size() == 0);
@@ -1457,10 +1481,12 @@ namespace {
 
 
 
+    template<class T>
     std::vector<ADB>
-    FullyImplicitBlackoilSolver::computeRelPerm(const SolutionState& state) const
+    FullyImplicitBlackoilSolver<T>::computeRelPerm(const SolutionState& state) const
     {
-        const int               nc   = grid_.number_of_cells;
+        using namespace Opm::AutoDiffGrid;
+        const int               nc   = numCells(grid_);
         const std::vector<int>& bpat = state.pressure.blockPattern();
 
         const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
@@ -1482,10 +1508,12 @@ namespace {
     }
 
 
+    template<class T>
     std::vector<ADB>
-    FullyImplicitBlackoilSolver::computePressures(const SolutionState& state) const
+    FullyImplicitBlackoilSolver<T>::computePressures(const SolutionState& state) const
     {
-        const int               nc   = grid_.number_of_cells;
+        using namespace Opm::AutoDiffGrid;
+        const int               nc   = numCells(grid_);
         const std::vector<int>& bpat = state.pressure.blockPattern();
 
         const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
@@ -1523,8 +1551,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
     {
@@ -1554,8 +1583,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,
@@ -1595,8 +1625,9 @@ namespace {
 
 
 
+    template<class T>
     double
-    FullyImplicitBlackoilSolver::residualNorm() const
+    FullyImplicitBlackoilSolver<T>::residualNorm() const
     {
         double globalNorm = 0;
         std::vector<ADB>::const_iterator quantityIt = residual_.mass_balance.begin();
@@ -1620,8 +1651,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   ,
@@ -1645,8 +1677,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   ,
@@ -1670,8 +1703,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   ,
@@ -1696,8 +1730,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);
@@ -1707,15 +1742,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);
@@ -1725,8 +1762,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);
@@ -1734,8 +1772,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()) {
@@ -1761,8 +1800,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()) {
@@ -1786,8 +1826,9 @@ namespace {
 
 
     /*
+    template<class T>
     void
-    FullyImplicitBlackoilSolver::
+    FullyImplicitBlackoilSolver<T>::
     classifyCondition(const SolutionState&        state,
                       std::vector<PhasePresence>& cond ) const
     {
@@ -1827,10 +1868,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,7 +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
@@ -39,49 +41,52 @@ namespace Opm
 
         /// Construct contained derived geological properties
         /// from grid and property information.
-        template <class Props>
+        template <class Props, class Grid>
-        DerivedGeology(const UnstructuredGrid& grid,
+        DerivedGeology(const Grid&             grid,
                        const Props&            props ,
                        const double*           grav = 0)
-            : pvol_ (grid.number_of_cells)
+            : pvol_ (Opm::AutoDiffGrid::numCells(grid))
-            , trans_(grid.number_of_faces)
+            , trans_(Opm::AutoDiffGrid::numFaces(grid))
-            , gpot_ (Vector::Zero(grid.cell_facepos[ grid.number_of_cells ], 1))
+            , gpot_ (Vector::Zero(Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1))
-            , z_(grid.number_of_cells)
+            , z_(Opm::AutoDiffGrid::numCells(grid))
         {
+            using namespace Opm::AutoDiffGrid;
             // Pore volume
-            const typename Vector::Index nc = grid.number_of_cells;
+            const typename Vector::Index nc = numCells(grid);
-            std::transform(grid.cell_volumes, grid.cell_volumes + nc,
+            std::transform(beginCellVolumes(grid), endCellVolumes(grid),
                            props.porosity(), pvol_.data(),
                            std::multiplies<double>());
 
             // Transmissibility
-            Vector htrans(grid.cell_facepos[nc]);
+            Vector htrans(numCellFaces(grid));
-            UnstructuredGrid* ug = const_cast<UnstructuredGrid*>(& grid);
+            Grid* ug = const_cast<Grid*>(& grid);
             tpfa_htrans_compute(ug, props.permeability(), htrans.data());
             tpfa_trans_compute (ug, htrans.data()     , trans_.data());
 
             // Compute z coordinates
             for (int c = 0; c<nc; ++c){
-                z_[c] = grid.cell_centroids[c*3 + 2];
+                z_[c] = cellCentroid(grid, c)[2];
             }
 
 
             // Gravity potential
             std::fill(gravity_, gravity_ + 3, 0.0);
             if (grav != 0) {
-                const typename Vector::Index nd = grid.dimensions;
+                const typename Vector::Index nd = dimensions(grid);
+                typedef typename ADCell2FacesTraits<Grid>::Type Cell2Faces;
+                Cell2Faces c2f=cell2Faces(grid);
 
                 for (typename Vector::Index c = 0; c < nc; ++c) {
-                    const double* const cc = & grid.cell_centroids[c*nd + 0];
+                    const double* const cc = cellCentroid(grid, c);
 
-                    const int* const p = grid.cell_facepos;
+                    typename Cell2Faces::row_type faces=c2f[c];
-                    for (int i = p[c]; i < p[c + 1]; ++i) {
+                    typedef typename Cell2Faces::row_type::iterator Iter;
-                        const int f = grid.cell_faces[i];
 
-                        const double* const fc = & grid.face_centroids[f*nd + 0];
+                    for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f) {
+                        const double* const fc = faceCentroid(grid, *f);
 
                         for (typename Vector::Index d = 0; d < nd; ++d) {
-                            gpot_[i] += grav[d] * (fc[d] - cc[d]);
+                            gpot_[f-faces.begin()] += grav[d] * (fc[d] - cc[d]);
                         }
                     }
                 }

opm/autodiff/GridHelpers.cpp

@@ -0,0 +1,278 @@
+/*
+  Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services.
+  Copyright 2014 Statoil AS
+
+  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 "config.h"
+
+#include <opm/autodiff/GridHelpers.hpp>
+namespace Opm
+{
+namespace AutoDiffGrid
+{
+
+// Interface functions using Unstructured grid
+/*    
+int numCells(const UnstructuredGrid& grid)
+{
+    return grid.number_of_cells;
+}
+
+int numFaces(const UnstructuredGrid& grid)
+{
+    return grid.number_of_faces;
+}
+int dimensions(const UnstructuredGrid& grid)
+{
+    return grid.dimensions;
+}
+*/
+Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
+faceCells(const UnstructuredGrid& grid)
+{
+    typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
+    return Eigen::Map<TwoColInt>(grid.face_cells, grid.number_of_faces, 2);
+}
+
+Eigen::Array<double, Eigen::Dynamic, 1>
+cellCentroidsZ(const UnstructuredGrid& grid)
+{
+    return Eigen::Map<Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >
+        (grid.cell_centroids, grid.number_of_cells, grid.dimensions).rightCols<1>();
+}
+
+const double*
+cellCentroid(const UnstructuredGrid& grid, int cell_index)
+{
+    return grid.cell_centroids+(cell_index*grid.dimensions);
+}
+
+const double* faceCentroid(const UnstructuredGrid& grid, int face_index)
+{
+    return grid.face_centroids+(face_index*grid.dimensions);
+}
+/*
+SparseTableView cell2Faces(const UnstructuredGrid& grid)
+{
+    return SparseTableView(grid.cell_faces, grid.cell_facepos, numCells(grid));
+}
+*/
+double cellVolume(const UnstructuredGrid& grid, int cell_index)
+{
+    return grid.cell_volumes[cell_index];
+}
+
+const double* beginCellVolumes(const UnstructuredGrid& grid)
+{
+    return grid.cell_volumes;
+}
+const double* endCellVolumes(const UnstructuredGrid& grid)
+{
+    return grid.cell_volumes+numCells(grid);
+}
+
+void extractInternalFaces(const UnstructuredGrid& grid,
+                          Eigen::Array<int, Eigen::Dynamic, 1>& internal_faces,
+                          Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>& nbi)
+{
+    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 = faceCells(grid);
+    // std::cout << "nb = \n" << nb << std::endl;
+    // Extracts the internal faces of the grid.
+    // These are stored in internal_faces.
+    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;
+    nbi.resize(num_internal, 2);
+    internal_faces.resize(num_internal);
+    int fi = 0;
+    int nf = numFaces(grid);
+
+    for (int f = 0; f < nf; ++f) {
+        if (ifaces[f]) {
+            internal_faces[fi] = f;
+            nbi.row(fi) = nb.row(f);
+            ++fi;
+        }
+    }
+}
+} // end namespace AutoDiffGrid
+
+#ifdef HAVE_DUNE_CORNERPOINT
+// Interface functions using CpGrid
+
+namespace UgGridHelpers
+{
+
+int numCells(const Dune::CpGrid& grid)
+{
+    return grid.numCells();
+}
+
+int numFaces(const  Dune::CpGrid& grid)
+{
+    return grid.numFaces();
+}
+
+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]);
+}
+
+const int*  globalCell(const Dune::CpGrid& grid)
+{
+    return &(grid.globalCell()[0]);
+}
+
+CellCentroidTraits<Dune::CpGrid>::IteratorType
+beginCellCentroids(const Dune::CpGrid& 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]);
+}
+
+double faceArea(const Dune::CpGrid& grid, int face_index)
+{
+    return grid.faceArea(face_index);
+}
+} // end namespace UgGridHelpers
+
+namespace AutoDiffGrid
+{
+
+
+Eigen::Array<double, Eigen::Dynamic, 1>
+cellCentroidsZ(const Dune::CpGrid& grid)
+{
+    // Create an Eigen array of appropriate size
+    int rows=numCells(grid);
+    Eigen::Array<double, Eigen::Dynamic, 1> array(rows);
+    // Fill it with the z coordinate of the cell centroids.
+    for (int i=0; i<rows; ++i)
+        array[i]=cellCentroid(grid, i)[2];
+    return array;
+}
+
+const double* cellCentroid(const Dune::CpGrid& grid, int cell_index)
+{
+    return &(grid.cellCentroid(cell_index)[0]);
+}
+
+const double* faceCentroid(const Dune::CpGrid& grid, int face_index)
+{
+    return &(grid.faceCentroid(face_index)[0]);
+}
+
+double cellVolume(const  Dune::CpGrid& grid, int cell_index)
+{
+    return grid.cellVolume(cell_index);
+}
+
+void extractInternalFaces(const Dune::CpGrid& grid,
+                          Eigen::Array<int, Eigen::Dynamic, 1>& internal_faces,
+                          Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>& nbi)
+{
+    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;
+    ADFaceCellTraits<Dune::CpGrid>::Type nb = faceCells(grid);
+    // std::cout << "nb = \n" << nb << std::endl;
+    // Extracts the internal faces of the grid.
+    // These are stored in internal_faces.
+    int nf=numFaces(grid);
+    int num_internal=0;
+    for(int f=0; f<nf; ++f)
+    {
+        if(grid.faceCell(f, 0)<0 || grid.faceCell(f, 1)<0)
+            continue;
+        ++num_internal;
+    }
+    // std::cout << num_internal << " internal faces." << std::endl;
+    nbi.resize(num_internal, 2);
+    internal_faces.resize(num_internal);
+    int fi = 0;
+
+    for (int f = 0; f < nf; ++f) {
+        if(grid.faceCell(f, 0)>=0 && grid.faceCell(f, 1)>=0) {
+            internal_faces[fi] = f;
+            nbi(fi,0) = grid.faceCell(f, 0);
+            nbi(fi,1) = grid.faceCell(f, 1);
+            ++fi;
+        }
+    }
+}
+
+
+CellVolumeIterator beginCellVolumes(const Dune::CpGrid& grid)
+{
+    return CellVolumeIterator(grid, 0);
+}
+
+CellVolumeIterator endCellVolumes(const Dune::CpGrid& grid)
+{
+    return CellVolumeIterator(grid, numCells(grid));
+}
+}       // end namespace AutoDiffGrid
+#endif  // HAVE_DUNE_CORNERPOINT
+}       // end namespace Opm

opm/autodiff/GridHelpers.hpp

@@ -0,0 +1,528 @@
+/*
+  Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services.
+  Copyright 2014 Statoil AS
+
+  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_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>
+#include <Eigen/Eigen>
+#include <Eigen/Sparse>
+
+#ifdef HAVE_DUNE_CORNERPOINT
+#include <dune/grid/CpGrid.hpp>
+#endif
+
+namespace Opm
+{
+
+namespace UgGridHelpers
+{
+
+} //end namespace UgGridHelpers
+namespace AutoDiffGrid
+{
+
+/// \brief Mapps a grid type to the corresponding face to cell mapping.
+///
+/// The value of the mapping is provided by the type Type.
+template<class T>
+struct ADFaceCellTraits
+{
+};
+
+/// \brief Get the z coordinates of the cell centroids of a grid.
+Eigen::Array<double, Eigen::Dynamic, 1>
+cellCentroidsZ(const UnstructuredGrid& 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 UnstructuredGrid& 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 UnstructuredGrid& grid, int face_index);
+
+/// \brief Mapping of the grid type to the type of the cell to faces mapping.
+template<class T>
+struct ADCell2FacesTraits
+    : public Opm::UgGridHelpers::Cell2FacesTraits<T>
+{
+};
+
+/// \brief Get the volume of a cell.
+/// \param grid The grid the cell belongs to.
+/// \param cell_index The index of the cell.
+double cellVolume(const UnstructuredGrid& grid, int cell_index);
+
+/// \brief The mapping of the grid type to type of the iterator over
+/// the cell volumes.
+///
+/// The value of the mapping is stored in nested type IteratorType
+/// \tparam T The type of the grid.
+template<class T>
+struct ADCellVolumesTraits
+{
+};
+
+template<>
+struct ADCellVolumesTraits<UnstructuredGrid>
+{
+    typedef const double* IteratorType;
+};
+
+/// \brief Get an iterator over the cell volumes of a grid positioned at the first cell.
+const double* beginCellVolumes(const UnstructuredGrid& grid);
+
+/// \brief Get an iterator over the cell volumes of a grid positioned after the last cell.
+const double* endCellVolumes(const UnstructuredGrid& 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 UnstructuredGrid& grid,
+                          Eigen::Array<int, Eigen::Dynamic, 1>& internal_faces,
+                          Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>& nbi);
+
+} // end namespace AutoDiffGrid
+} // end namespace Opm
+
+#ifdef HAVE_DUNE_CORNERPOINT
+
+#include <dune/common/iteratorfacades.hh>
+namespace Opm
+{
+
+namespace AutoDiffGrid
+{
+/// \brief A proxy class representing a row of FaceCellsContainer.
+class FaceCellsProxy
+{    
+public:
+    /// \brief Constructor.
+    /// \param grid The grid whose face to cell mapping we represent.
+    /// \param cell_index The index of the cell we repesent.
+    FaceCellsProxy(const Dune::CpGrid* grid, int cell_index)
+        : grid_(grid), cell_index_(cell_index)
+    {}
+    /// \brief Get the index of the cell associated with a local_index.
+    int operator[](int local_index)
+    {
+        return grid_->faceCell(cell_index_, local_index);
+    }
+private:
+    const Dune::CpGrid* grid_;
+    int cell_index_;
+};
+
+/// \brief A class representing the face to cells mapping similar to the
+/// way done in UnstructuredGrid.
+class FaceCellsContainerProxy
+{
+public:
+    typedef FaceCellsProxy row_type;
+    
+    /// \brief Constructor.
+    /// \param grid The grid whose information we represent.
+    FaceCellsContainerProxy(const Dune::CpGrid* grid)
+        : grid_(grid)
+    {}
+    /// \brief Get the mapping for a cell.
+    /// \param cell_index The index of the cell.
+    FaceCellsProxy operator[](int cell_index) const
+    {
+        return FaceCellsProxy(grid_, cell_index);
+    }
+    /// \brief Get a face associated with a cell.
+    /// \param cell_index The index of the cell.
+    /// \param local_index The local index of the cell, either 0 or 1.
+    /// \param The index of the face or -1 if it is not present because of
+    /// a boundary.
+    int operator()(int cell_index, int local_index) const
+    {
+        return grid_->faceCell(cell_index, local_index);
+    }
+private:
+    const Dune::CpGrid* grid_;
+};
+
+class Cell2FacesRow
+{
+public:
+    class iterator
+        : public Dune::RandomAccessIteratorFacade<iterator,int, int, int>
+    {
+    public:
+        iterator(const Dune::cpgrid::OrientedEntityTable<0,1>::row_type* row,
+                 int index)
+            : row_(row), index_(index)
+        {}
+
+        void increment()
+        {
+            ++index_;
+        }
+        void decrement()
+        {
+            --index_;
+        }
+        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_;
+    };
+    
+    typedef iterator const_iterator;
+
+    Cell2FacesRow(const Dune::cpgrid::OrientedEntityTable<0,1>::row_type row)
+        : row_(row)
+    {}
+
+    const_iterator begin() const
+    {
+        return const_iterator(&row_, 0);
+    }
+
+    const_iterator end() const
+    {
+        return const_iterator(&row_, row_.size());
+    }
+    
+private:
+    const Dune::cpgrid::OrientedEntityTable<0,1>::row_type row_;
+};
+
+class Cell2FacesContainer
+{
+public:
+    typedef  Cell2FacesRow row_type;
+    
+    Cell2FacesContainer(const Dune::CpGrid* grid)
+        : grid_(grid)
+    {};
+    
+    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_, n);
+    }
+    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);
+
+double faceArea(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>
+{
+    typedef Cell2FacesContainer Type;
+};
+
+/// \brief Get the volume of a cell.
+/// \param grid The grid the cell belongs to.
+/// \param cell_index The index of the cell.
+double cellVolume(const  Dune::CpGrid& grid, int cell_index);
+
+/// \brief An iterator over the cell volumes.
+class CellVolumeIterator
+    : public Dune::RandomAccessIteratorFacade<CellVolumeIterator, double, double, int>
+{
+public:
+    /// \brief Creates an iterator.
+    /// \param grid The grid the iterator belongs to.
+    /// \param cell_index The position of the iterator.
+    CellVolumeIterator(const  Dune::CpGrid& grid, int cell_index)
+        : grid_(&grid), cell_index_(cell_index)
+    {}
+
+    double dereference() const
+    {
+        return grid_->cellVolume(cell_index_);
+    }
+    void increment()
+    {
+        ++cell_index_;
+    }
+    double elementAt(int n) const
+    {
+        return grid_->cellVolume(n);
+    }
+    void advance(int n)
+    {
+        cell_index_+=n;
+    }
+    void decrement()
+    {
+        --cell_index_;
+    }
+    int distanceTo(const CellVolumeIterator& o) const
+    {
+        return o.cell_index_-cell_index_;
+    }
+    bool equals(const CellVolumeIterator& o) const
+    {
+        return o.grid_==grid_ && o.cell_index_==cell_index_;
+    }
+    
+private:
+    const Dune::CpGrid* grid_;
+    int cell_index_;
+};
+
+template<>
+struct ADCellVolumesTraits<Dune::CpGrid>
+{
+    typedef CellVolumeIterator IteratorType;
+};
+
+/// \brief Get an iterator over the cell volumes of a grid positioned at the first cell.
+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,9 +16,11 @@
   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>
+#include <opm/autodiff/GridHelpers.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/simulator/WellState.hpp>
@@ -54,15 +56,24 @@ namespace {
                     const HelperOps&        ops ,
                     const GeoProps&         geo )
     {
-        const int nc = grid.number_of_cells;
+        using namespace Opm::AutoDiffGrid;
+        const int nc = numCells(grid);
 
-        std::vector<int> f2hf(2 * grid.number_of_faces, -1);
+        std::vector<int> f2hf(2 * numFaces(grid), -1);
-        for (int c = 0, i = 0; c < nc; ++c) {
+        Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
-            for (; i < grid.cell_facepos[c + 1]; ++i) {
+            face_cells;
-                const int f = grid.cell_faces[ i ];
+        
-                const int p = 0 + (grid.face_cells[2*f + 0] != c);
+        typedef typename Opm::UgGridHelpers::Cell2FacesTraits<UnstructuredGrid>::Type
-
+            Cell2Faces;
-                f2hf[2*f + p] = i;
+        Cell2Faces c2f=cell2Faces(grid);
+        for (int c = 0; c < nc; ++c) {
+            typename Cell2Faces::row_type 
+                cell_faces = c2f[c];
+            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);
+                f2hf[2*(*f) + p] = f-c2f[0].begin();
             }
         }
 
@@ -78,8 +89,10 @@ namespace {
         std::vector<Tri> grav;  grav.reserve(2 * ni);
         for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
             const int f  = ops.internal_faces[ i ];
-            const int c1 = grid.face_cells[2*f + 0];
+            Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
-            const int c2 = grid.face_cells[2*f + 1];
+                face_cells=faceCells(grid);
+            const int c1 = face_cells(f,0);
+            const int c2 = face_cells(f,1);
 
             assert ((c1 >= 0) && (c2 >= 0));
 
@@ -98,9 +111,10 @@ namespace {
 
     V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav)
     {
+        using namespace Opm::AutoDiffGrid;
         const int nw = wells.number_of_wells;
         const int nperf = wells.well_connpos[nw];
-        const int dim = grid.dimensions;
+        const int dim = dimensions(grid);
         V wdp = V::Zero(nperf,1);
         assert(wdp.size() == rho.size());
 
@@ -157,7 +171,8 @@ namespace {
                        BlackoilState& state,
                        WellState& well_state)
     {
-        const int nc = grid_.number_of_cells;
+        using namespace Opm::AutoDiffGrid;
+        const int nc = numCells(grid_);
         const int np = state.numPhases();
 
         well_flow_residual_.resize(np, ADB::null());
@@ -226,15 +241,16 @@ namespace {
                                      const BlackoilState& state,
                                      const WellState& well_state)
     {
+        using namespace Opm::AutoDiffGrid;
         // Suppress warnings about "unused parameters".
         static_cast<void>(dt);
         static_cast<void>(well_state);
 
-        const int nc = grid_.number_of_cells;
+        const int nc = numCells(grid_);
         const int np = state.numPhases();
         const int nw = wells_.number_of_wells;
         const int nperf = wells_.well_connpos[nw];
-        const int dim = grid_.dimensions;
+        const int dim = dimensions(grid_);
 
         const std::vector<int> cells = buildAllCells(nc);
 
@@ -284,9 +300,9 @@ namespace {
                           const BlackoilState& state,
                           const WellState& well_state)
     {
-
+        using namespace Opm::AutoDiffGrid;
         const V& pv = geo_.poreVolume();
-        const int nc = grid_.number_of_cells;
+        const int nc = numCells(grid_);              ;
         const int np = state.numPhases();
         const int nw = wells_.number_of_wells;
         const int nperf = wells_.well_connpos[nw];
@@ -415,7 +431,8 @@ namespace {
     ImpesTPFAAD::solveJacobianSystem(BlackoilState& state,
                                      WellState& well_state) const
     {
-        const int nc = grid_.number_of_cells;
+        using namespace Opm::AutoDiffGrid;
+        const int nc = numCells(grid_);
         const int nw = wells_.number_of_wells;
         // const int np = state.numPhases();
 
@@ -458,9 +475,10 @@ namespace {
     ImpesTPFAAD::computeFluxes(BlackoilState& state,
                                WellState& well_state) const
     {
+        using namespace Opm::AutoDiffGrid;
         // This method computes state.faceflux(),
         // well_state.perfRates() and well_state.perfPress().
-        const int nc = grid_.number_of_cells;
+        const int nc = numCells(grid_);
         const int np = state.numPhases();
         const int nw = wells_.number_of_wells;
         const int nperf = wells_.well_connpos[nw];
@@ -515,8 +533,8 @@ namespace {
             flux += face_mob * head;
         }
 
-        V all_flux = superset(flux, ops_.internal_faces, grid_.number_of_faces);
+        V all_flux = superset(flux, ops_.internal_faces, numFaces(grid_));
-        std::copy(all_flux.data(), all_flux.data() + grid_.number_of_faces, state.faceflux().begin());
+        std::copy(all_flux.data(), all_flux.data() + numFaces(grid_), state.faceflux().begin());
 
         perf_flux = -perf_flux; // well_state.perfRates() assumed to be inflows.
         std::copy(perf_flux.data(), perf_flux.data() + nperf, well_state.perfRates().begin());

opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -40,9 +40,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
@@ -66,7 +69,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,
@@ -92,4 +95,5 @@ namespace Opm
 
 } // namespace Opm
 
+#include "SimulatorFullyImplicitBlackoil_impl.hpp"
 #endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp

@@ -0,0 +1,196 @@
+#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>
+
+#ifdef HAVE_DUNE_CORNERPOINT
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#endif
+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)
+    {
+        // Write data in VTK format.
+        std::ostringstream vtkfilename;
+        std::ostringstream vtkpath;
+        vtkpath << output_dir << "/vtk_files";
+        vtkpath << "/output-" << std::setw(3) << std::setfill('0') << step;
+        boost::filesystem::path fpath(vtkpath.str());
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        vtkfilename << "output-" << std::setw(3) << std::setfill('0') << step;
+        Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafView(), Dune::VTK::nonconforming);
+        writer.addCellData(state.saturation(), "saturation", state.numPhases());
+        writer.addCellData(state.pressure(), "pressure", 1);
+        
+        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);
+        writer.addCellData(cell_velocity, "velocity", Dune::CpGrid::dimension);
+        writer.pwrite(vtkfilename.str(), vtkpath.str(), std::string("."), Dune::VTK::ascii);
+    }
+#endif
+
+}

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -0,0 +1,91 @@
+#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 <opm/core/utility/DataMap.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+
+#include <opm/autodiff/GridHelpers.hpp>
+
+#include <string>
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+
+#include <boost/filesystem.hpp>
+
+#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);
+#endif
+
+    template<class Grid>
+    void outputStateMatlab(const Grid& 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"));
+        }
+    }
+
+
+}
+#endif

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -17,11 +17,7 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-
+#include  <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
-#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>
@@ -30,7 +26,6 @@
 #include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
-
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
 #include <opm/core/pressure/flow_bc.h>
@@ -61,12 +56,12 @@
 
 namespace Opm
 {
-
+    template<class T>
-    class SimulatorFullyImplicitBlackoil::Impl
+    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,
@@ -89,7 +84,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_;
@@ -97,7 +92,7 @@ namespace Opm
         const double* gravity_;
         // Solvers
         DerivedGeology geo_;
-        FullyImplicitBlackoilSolver solver_;
+        FullyImplicitBlackoilSolver<Grid> solver_;
         // Misc. data
         std::vector<int> allcells_;
     };
@@ -105,8 +100,9 @@ namespace Opm
 
 
 
-    SimulatorFullyImplicitBlackoil::SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
+    template<class T>
-                                                                   const UnstructuredGrid& grid,
+    SimulatorFullyImplicitBlackoil<T>::SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
+                                                                   const Grid& grid,
                                                                    BlackoilPropsAdInterface& props,
                                                                    const RockCompressibility* rock_comp_props,
                                                                    WellsManager& wells_manager,
@@ -121,7 +117,8 @@ namespace Opm
 
 
 
-    SimulatorReport SimulatorFullyImplicitBlackoil::run(SimulatorTimer& timer,
+    template<class T>
+    SimulatorReport SimulatorFullyImplicitBlackoil<T>::run(SimulatorTimer& timer,
                                                         BlackoilState& state,
                                                         WellStateFullyImplicitBlackoil& well_state)
     {
@@ -130,70 +127,6 @@ 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(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(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::WellStateFullyImplicitBlackoil& well_state,
                                   const int step,
                                   const std::string& output_dir)
@@ -252,8 +185,9 @@ namespace Opm
 
 
     // \TODO: Treat bcs.
-    SimulatorFullyImplicitBlackoil::Impl::Impl(const parameter::ParameterGroup& param,
+    template<class T>
-                                               const UnstructuredGrid& grid,
+    SimulatorFullyImplicitBlackoil<T>::Impl::Impl(const parameter::ParameterGroup& param,
+                                               const Grid& grid,
                                                BlackoilPropsAdInterface& props,
                                                const RockCompressibility* rock_comp_props,
                                                WellsManager& wells_manager,
@@ -293,26 +227,24 @@ namespace Opm
         max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
 
         // Misc init.
-        const int num_cells = grid.number_of_cells;
+        const int num_cells = AutoDiffGrid::numCells(grid);
         allcells_.resize(num_cells);
         for (int cell = 0; cell < num_cells; ++cell) {
             allcells_[cell] = cell;
         }
     }
 
-
+    template<class T>
-
+    SimulatorReport SimulatorFullyImplicitBlackoil<T>::Impl::run(SimulatorTimer& timer,
-
-    SimulatorReport SimulatorFullyImplicitBlackoil::Impl::run(SimulatorTimer& timer,
                                                               BlackoilState& state,
                                                               WellStateFullyImplicitBlackoil& well_state)
     {
         // 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;
@@ -392,7 +324,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

opm/autodiff/SimulatorIncompTwophaseAd.cpp

@@ -23,6 +23,7 @@
 #endif // HAVE_CONFIG_H
 
 #include <opm/autodiff/SimulatorIncompTwophaseAd.hpp>
+#include <opm/autodiff/GridHelpers.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 
@@ -400,7 +401,7 @@ namespace Opm
         num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
 
         // Misc init.
-        const int num_cells = grid.number_of_cells;
+        const int num_cells = Opm::AutoDiffGrid::numCells(grid);
         allcells_.resize(num_cells);
         for (int cell = 0; cell < num_cells; ++cell) {
             allcells_[cell] = cell;
@@ -498,11 +499,13 @@ namespace Opm
                     double av_prev_press = 0.0;
                     double av_press = 0.0;
                     double tot_vol = 0.0;
-                    const int num_cells = grid_.number_of_cells;
+                    const int num_cells = Opm::AutoDiffGrid::numCells(grid_);
                     for (int cell = 0; cell < num_cells; ++cell) {
-                        av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
+                        av_prev_press += initial_pressure[cell]*
-                        av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
+                            Opm::AutoDiffGrid::cellVolume(grid_, cell);
-                        tot_vol       += grid_.cell_volumes[cell];
+                        av_press      += state.pressure()[cell]*
+                            Opm::AutoDiffGrid::cellVolume(grid_, cell);
+                        tot_vol       += Opm::AutoDiffGrid::cellVolume(grid_, cell);
                     }
                     // Renormalization constant
                     const double ren_const = (av_prev_press - av_press)/tot_vol;

opm/autodiff/TransportSolverTwophaseAd.cpp

@@ -51,21 +51,22 @@ namespace Opm
           tol_(param.getDefault("nl_tolerance", 1e-9)),
           maxit_(param.getDefault("nl_maxiter", 30))
     {
-        const int nc = grid_.number_of_cells;
+        using namespace Opm::AutoDiffGrid;
+        const int nc = numCells(grid_);
         allcells_.resize(nc);
         for (int i = 0; i < nc; ++i) {
             allcells_[i] = i;
         }
-        if (gravity && gravity[grid_.dimensions - 1] != 0.0) {
+        if (gravity && gravity[dimensions(grid_) - 1] != 0.0) {
-            gravity_ = gravity[grid_.dimensions - 1];
+            gravity_ = gravity[dimensions(grid_) - 1];
-            for (int dd = 0; dd < grid_.dimensions - 1; ++dd) {
+            for (int dd = 0; dd < dimensions(grid_) - 1; ++dd) {
                 if (gravity[dd] != 0.0) {
                     OPM_THROW(std::runtime_error, "TransportSolverTwophaseAd: can only handle gravity aligned with last dimension");
                 }
             }
             V htrans(grid.cell_facepos[grid.number_of_cells]);
             tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid), props.permeability(), htrans.data());
-            V trans(grid_.number_of_faces);
+            V trans(numFaces(grid_));
             tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid), htrans.data(), trans.data());
             transi_ = subset(trans, ops_.internal_faces);
         }
@@ -161,16 +162,17 @@ namespace Opm
                                           const double dt,
                                           TwophaseState& state)
     {
+        using namespace Opm::AutoDiffGrid;
         typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
         typedef Eigen::Map<const V> Vec;
-        const int nc = grid_.number_of_cells;
+        const int nc = numCells(grid_);
         const TwoCol s0 = Eigen::Map<const TwoCol>(state.saturation().data(), nc, 2);
         double res_norm = 1e100;
         const V sw0 = s0.leftCols<1>();
         // sw1 is the object that will be changed every Newton iteration.
         // V sw1 = sw0;
         V sw1 = 0.5*V::Ones(nc,1);
-        const V dflux_all = Vec(state.faceflux().data(), grid_.number_of_faces, 1);
+        const V dflux_all = Vec(state.faceflux().data(), numFaces(grid_), 1);
         const int num_internal = ops_.internal_faces.size();
         V dflux = subset(dflux_all, ops_.internal_faces);
 
@@ -184,7 +186,7 @@ namespace Opm
         const V p1 = Vec(state.pressure().data(), nc, 1);
         const V ndp = (ops_.ngrad * p1.matrix()).array();
         typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> DynArr;
-        const V z = Eigen::Map<DynArr>(grid_.cell_centroids, nc, grid_.dimensions).rightCols<1>();
+        const V z = cellCentroidsZ(grid_);
         const V ndz = (ops_.ngrad * z.matrix()).array();
         assert(num_internal == ndp.size());
         const double* density = props_.density();