Merge pull request #544 from osae/equil-init

Equil init

examples/compute_initial_state.cpp

@@ -0,0 +1,110 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/simulator/initStateEquil.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
+#include <boost/filesystem.hpp>
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "----------------------------------------------------------------" << std::endl;
+        }
+    }
+
+    void outputData(const std::string& output_dir,
+                    const std::string& name,
+                    const std::vector<double>& data)
+    {
+        std::ostringstream fname;
+        fname << output_dir << "/" << name;
+        boost::filesystem::path fpath = fname.str();
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        fname << "/" << "initial.txt";
+        std::ofstream file(fname.str().c_str());
+        if (!file) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
+        }
+        std::copy(data.begin(), data.end(), std::ostream_iterator<double>(file, "\n"));
+    }
+
+
+
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    using namespace Opm;
+
+    // Setup.
+    parameter::ParameterGroup param(argc, argv, false);
+    std::cout << "---------------    Reading parameters     ---------------" << std::endl;
+    const std::string deck_filename = param.get<std::string>("deck_filename");
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile(deck_filename);
+    const double grav = param.getDefault("gravity", unit::gravity);
+    //EclipseGridParser deck(deck_filename);
+    GridManager gm(deck);
+    const UnstructuredGrid& grid = *gm.c_grid();
+    BlackoilPropertiesFromDeck props(deck, grid, param);
+    warnIfUnusedParams(param);
+
+    // Initialisation.
+    BlackoilState state;
+    initStateEquil(grid, props, deck, grav, state);
+
+    // Output.
+    const std::string output_dir = param.getDefault<std::string>("output_dir", "output");
+    outputData(output_dir, "pressure", state.pressure());
+    outputData(output_dir, "saturation", state.saturation());
+    outputData(output_dir, "rs", state.gasoilratio());
+    outputData(output_dir, "rv", state.rv());
+}
+catch (const std::exception& e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}

opm/core/simulator/EquilibrationHelpers.hpp

@@ -0,0 +1,829 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
+#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
+
+#include <opm/core/props/BlackoilPropertiesInterface.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/core/utility/linearInterpolation.hpp>
+#include <opm/core/utility/RegionMapping.hpp>
+#include <opm/core/utility/RootFinders.hpp>
+
+#include <memory>
+
+
+/*
+---- synopsis of EquilibrationHelpers.hpp ----
+
+namespace Opm
+{
+    namespace Equil {
+
+        template <class Props>
+        class DensityCalculator;
+
+        template <>
+        class DensityCalculator< BlackoilPropertiesInterface >;
+
+        namespace Miscibility {
+            class RsFunction;
+            class NoMixing;
+            class RsVD;
+            class RsSatAtContact;
+        }
+
+        struct EquilRecord;
+
+        template <class DensCalc>
+        class EquilReg;
+
+
+        struct PcEq;
+
+        inline double satFromPc(const BlackoilPropertiesInterface& props,
+                                const int phase,
+                                const int cell,
+                                const double target_pc,
+                                const bool increasing = false);
+        struct PcEqSum
+        inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
+                                      const int phase1,
+                                      const int phase2,
+                                      const int cell,
+                                      const double target_pc);
+    } // namespace Equil
+} // namespace Opm
+
+---- end of synopsis of EquilibrationHelpers.hpp ----
+*/
+
+
+namespace Opm
+{
+    /**
+     * Types and routines that collectively implement a basic
+     * ECLIPSE-style equilibration-based initialisation scheme.
+     *
+     * This namespace is intentionally nested to avoid name clashes
+     * with other parts of OPM.
+     */
+    namespace Equil {
+
+
+        template <class Props>
+        class DensityCalculator;
+
+        /**
+         * Facility for calculating phase densities based on the
+         * BlackoilPropertiesInterface.
+         *
+         * Implements the crucial <CODE>operator()(p,svol)</CODE>
+         * function that is expected by class EquilReg.
+         */
+        template <>
+        class DensityCalculator< BlackoilPropertiesInterface > {
+        public:
+            /**
+             * Constructor.
+             *
+             * \param[in] props Implementation of the
+             * BlackoilPropertiesInterface.
+             *
+             * \param[in] c Single cell used as a representative cell
+             * in a PVT region.
+             */
+            DensityCalculator(const BlackoilPropertiesInterface& props,
+                              const int                          c)
+                : props_(props)
+                , c_(1, c)
+            {
+            }
+
+            /**
+             * Compute phase densities of all phases at phase point
+             * given by (pressure, surface volume) tuple.
+             *
+             * \param[in] p Fluid pressure.
+             *
+             * \param[in] z Surface volumes of all phases.
+             *
+             * \return Phase densities at phase point.
+             */
+            std::vector<double>
+            operator()(const double               p,
+                       const std::vector<double>& z) const
+            {
+                const int np = props_.numPhases();
+                std::vector<double> A(np * np, 0);
+
+                assert (z.size() == std::vector<double>::size_type(np));
+
+                double* dAdp = 0;
+                props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
+
+                std::vector<double> rho(np, 0.0);
+                props_.density(1, &A[0], &rho[0]);
+
+                return rho;
+            }
+
+        private:
+            const BlackoilPropertiesInterface& props_;
+            const std::vector<int>             c_;
+        };
+
+
+        /**
+         * Types and routines relating to phase mixing in
+         * equilibration calculations.
+         */
+        namespace Miscibility {
+
+            /**
+             * Base class for phase mixing functions.
+             */
+            class RsFunction
+            {
+            public:
+                /**
+                 * Function call operator.
+                 *
+                 * \param[in] depth Depth at which to calculate RS
+                 * value.
+                 *
+                 * \param[in] press Pressure at which to calculate RS
+                 * value.
+                 *
+                 * \return Dissolved gas-oil ratio (RS) at depth @c
+                 * depth and pressure @c press.
+                 */
+                virtual double operator()(const double depth,
+                                          const double press,
+                                          const double sat = 0.0) const = 0;
+            };
+
+
+            /**
+             * Type that implements "no phase mixing" policy.
+             */
+            class NoMixing : public RsFunction {
+            public:
+                /**
+                 * Function call.
+                 *
+                 * \param[in] depth Depth at which to calculate RS
+                 * value.
+                 *
+                 * \param[in] press Pressure at which to calculate RS
+                 * value.
+                 *
+                 * \return Dissolved gas-oil ratio (RS) at depth @c
+                 * depth and pressure @c press.  In "no mixing
+                 * policy", this is identically zero.
+                 */
+                double
+                operator()(const double /* depth */,
+                           const double /* press */,
+                           const double sat = 0.0) const
+                {
+                    return 0.0;
+                }
+            };
+
+
+            /**
+             * Type that implements "dissolved gas-oil ratio"
+             * tabulated as a function of depth policy.  Data
+             * typically taken from keyword 'RSVD'.
+             */
+            class RsVD : public RsFunction {
+            public:
+                /**
+                 * Constructor.
+                 *
+                 * \param[in] props      property object
+                 * \param[in] cell       any cell in the pvt region
+                 * \param[in] depth Depth nodes.
+                 * \param[in] rs Dissolved gas-oil ratio at @c depth.
+                 */
+                RsVD(const BlackoilPropertiesInterface& props,
+                     const int cell,
+                     const std::vector<double>& depth,
+                     const std::vector<double>& rs)
+                    : props_(props) 
+                    , cell_(cell)
+                    , depth_(depth)
+                    , rs_(rs)
+                {
+                    auto pu = props_.phaseUsage();
+                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
+                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                }
+
+                /**
+                 * Function call.
+                 *
+                 * \param[in] depth Depth at which to calculate RS
+                 * value.
+                 *
+                 * \param[in] press Pressure at which to calculate RS
+                 * value.
+                 *
+                 * \return Dissolved gas-oil ratio (RS) at depth @c
+                 * depth and pressure @c press.
+                 */
+                double
+                operator()(const double depth,
+                           const double press,
+                           const double sat_gas = 0.0) const
+                {
+                    if (sat_gas > 0.0) {
+                        return satRs(press);
+                    } else {
+                        return std::min(satRs(press), linearInterpolation(depth_, rs_, depth));
+                    }
+                }
+
+            private:
+                const BlackoilPropertiesInterface& props_;
+                const int cell_;
+                std::vector<double> depth_; /**< Depth nodes */
+                std::vector<double> rs_;    /**< Dissolved gas-oil ratio */
+                double z_[BlackoilPhases::MaxNumPhases];
+                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+                double satRs(const double press) const
+                {
+                    props_.matrix(1, &press, z_, &cell_, A_, 0);
+                    // Rs/Bo is in the gas row and oil column of A_.
+                    // 1/Bo is in the oil row and column.
+                    // Recall also that it is stored in column-major order.
+                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+                    const int np = props_.numPhases();
+                    return A_[np*opos + gpos] / A_[np*opos + opos];
+                }
+            };
+
+
+            /**
+             * Type that implements "vaporized oil-gas ratio"
+             * tabulated as a function of depth policy.  Data
+             * typically taken from keyword 'RVVD'.
+             */
+            class RvVD : public RsFunction {
+            public:
+                /**
+                 * Constructor.
+                 *
+                 * \param[in] props      property object
+                 * \param[in] cell       any cell in the pvt region
+                 * \param[in] depth Depth nodes.
+                 * \param[in] rv Dissolved gas-oil ratio at @c depth.
+                 */
+                RvVD(const BlackoilPropertiesInterface& props,
+                     const int cell,
+                     const std::vector<double>& depth,
+                     const std::vector<double>& rv)
+                    : props_(props) 
+                    , cell_(cell)
+                    , depth_(depth)
+                    , rv_(rv)
+                {
+                    auto pu = props_.phaseUsage();
+                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
+                }
+
+                /**
+                 * Function call.
+                 *
+                 * \param[in] depth Depth at which to calculate RV
+                 * value.
+                 *
+                 * \param[in] press Pressure at which to calculate RV
+                 * value.
+                 *
+                 * \return Vaporized oil-gas ratio (RV) at depth @c
+                 * depth and pressure @c press.
+                 */
+                double
+                operator()(const double depth,
+                           const double press,
+                           const double sat_oil = 0.0 ) const
+                {
+                    if (sat_oil > 0.0) {
+                        return satRv(press);
+                    } else {
+                        return std::min(satRv(press), linearInterpolation(depth_, rv_, depth));
+                    }
+                }
+
+            private:
+                const BlackoilPropertiesInterface& props_;
+                const int cell_;
+                std::vector<double> depth_; /**< Depth nodes */
+                std::vector<double> rv_;    /**< Vaporized oil-gas ratio */
+                double z_[BlackoilPhases::MaxNumPhases];
+                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+                double satRv(const double press) const
+                {
+                    props_.matrix(1, &press, z_, &cell_, A_, 0);
+                    // Rv/Bg is in the oil row and gas column of A_.
+                    // 1/Bg is in the gas row and column.
+                    // Recall also that it is stored in column-major order.
+                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+                    const int np = props_.numPhases();
+                    return A_[np*gpos + opos] / A_[np*gpos + gpos];
+                }
+            };
+
+
+            /**
+             * Class that implements "dissolved gas-oil ratio" (Rs)
+             * as function of depth and pressure as follows:
+             *
+             *   1. The Rs at the gas-oil contact is equal to the
+             *      saturated Rs value, Rs_sat_contact.
+             *
+             *   2. The Rs elsewhere is equal to Rs_sat_contact, but
+             *      constrained to the saturated value as given by the
+             *      local pressure.
+             *
+             * This should yield Rs-values that are constant below the
+             * contact, and decreasing above the contact.
+             */
+            class RsSatAtContact : public RsFunction {
+            public:
+                /**
+                 * Constructor.
+                 *
+                 * \param[in] props      property object
+                 * \param[in] cell       any cell in the pvt region
+                 * \param[in] p_contact  oil pressure at the contact
+                 */
+                RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
+                    : props_(props), cell_(cell)
+                {
+                    auto pu = props_.phaseUsage();
+                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
+                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                    rs_sat_contact_ = satRs(p_contact);
+                }
+
+                /**
+                 * Function call.
+                 *
+                 * \param[in] depth Depth at which to calculate RS
+                 * value.
+                 *
+                 * \param[in] press Pressure at which to calculate RS
+                 * value.
+                 *
+                 * \return Dissolved gas-oil ratio (RS) at depth @c
+                 * depth and pressure @c press.
+                 */
+                double
+                operator()(const double /* depth */,
+                           const double press,
+                           const double sat_gas = 0.0) const
+                {                     
+                    if (sat_gas > 0.0) {
+                        return satRs(press);
+                    } else {
+                        return std::min(satRs(press), rs_sat_contact_);
+                    }
+                }
+
+            private:
+                const BlackoilPropertiesInterface& props_;
+                const int cell_;
+                double z_[BlackoilPhases::MaxNumPhases];
+                double rs_sat_contact_;
+                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+                double satRs(const double press) const
+                {
+                    props_.matrix(1, &press, z_, &cell_, A_, 0);
+                    // Rs/Bo is in the gas row and oil column of A_.
+                    // 1/Bo is in the oil row and column.
+                    // Recall also that it is stored in column-major order.
+                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+                    const int np = props_.numPhases();
+                    return A_[np*opos + gpos] / A_[np*opos + opos];
+                }
+            };
+
+
+            /**
+             * Class that implements "vaporized oil-gas ratio" (Rv)
+             * as function of depth and pressure as follows:
+             *
+             *   1. The Rv at the gas-oil contact is equal to the
+             *      saturated Rv value, Rv_sat_contact.
+             *
+             *   2. The Rv elsewhere is equal to Rv_sat_contact, but
+             *      constrained to the saturated value as given by the
+             *      local pressure.
+             *
+             * This should yield Rv-values that are constant below the
+             * contact, and decreasing above the contact.
+             */
+            class RvSatAtContact : public RsFunction {
+            public:
+                /**
+                 * Constructor.
+                 *
+                 * \param[in] props      property object
+                 * \param[in] cell       any cell in the pvt region
+                 * \param[in] p_contact  oil pressure at the contact
+                 */
+                RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
+                    : props_(props), cell_(cell)
+                {
+                    auto pu = props_.phaseUsage();
+                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
+                    rv_sat_contact_ = satRv(p_contact);
+                }
+
+                /**
+                 * Function call.
+                 *
+                 * \param[in] depth Depth at which to calculate RV
+                 * value.
+                 *
+                 * \param[in] press Pressure at which to calculate RV
+                 * value.
+                 *
+                 * \return Dissolved oil-gas ratio (RV) at depth @c
+                 * depth and pressure @c press.
+                 */
+                double
+                operator()(const double /*depth*/,
+                           const double press,
+                           const double sat_oil = 0.0) const
+                {
+                    if (sat_oil > 0.0) {
+                        return satRv(press);
+                    } else {
+                        return std::min(satRv(press), rv_sat_contact_);
+                    }
+                }
+
+            private:
+                const BlackoilPropertiesInterface& props_;
+                const int cell_;
+                double z_[BlackoilPhases::MaxNumPhases];
+                double rv_sat_contact_;
+                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+                double satRv(const double press) const
+                {
+                    props_.matrix(1, &press, z_, &cell_, A_, 0);
+                    // Rv/Bg is in the oil row and gas column of A_.
+                    // 1/Bg is in the gas row and column.
+                    // Recall also that it is stored in column-major order.
+                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+                    const int np = props_.numPhases();
+                    return A_[np*gpos + opos] / A_[np*gpos + gpos];
+                }
+            };
+
+        } // namespace Miscibility
+
+
+
+        /**
+         * Equilibration record.
+         *
+         * Layout and contents inspired by first six items of
+         * ECLIPSE's 'EQUIL' records.  This is the minimum amount of
+         * input data needed to define phase pressures in an
+         * equilibration region.
+         *
+         * Data consists of three pairs of depth and pressure values:
+         *   1. main
+         *     - @c depth Main datum depth.
+         *     - @c press Pressure at datum depth.
+         *
+         *   2. woc
+         *     - @c depth Depth of water-oil contact
+         *     - @c press water-oil capillary pressure at water-oil contact.
+         *       Capillary pressure defined as "P_oil - P_water".
+         *
+         *   3. goc
+         *     - @c depth Depth of gas-oil contact
+         *     - @c press Gas-oil capillary pressure at gas-oil contact.
+         *       Capillary pressure defined as "P_gas - P_oil".
+         */
+        struct EquilRecord {
+            struct {
+                double depth;
+                double press;
+            } main, woc, goc;
+            int live_oil_table_index;
+            int wet_gas_table_index;
+            int N;
+        };
+
+        /**
+         * Aggregate information base of an equilibration region.
+         *
+         * Provides inquiry methods for retrieving depths of contacs
+         * and pressure values as well as a means of calculating fluid
+         * densities, dissolved gas-oil ratio and vapourised oil-gas
+         * ratios.
+         *
+         * \tparam DensCalc Type that provides access to a phase
+         * density calculation facility.  Must implement an operator()
+         * declared as
+         * <CODE>
+         * std::vector<double>
+         * operator()(const double               press,
+         *            const std::vector<double>& svol )
+         * </CODE>
+         * that calculates the phase densities of all phases in @c
+         * svol at fluid pressure @c press.
+         */
+        template <class DensCalc>
+        class EquilReg {
+        public:
+            /**
+             * Constructor.
+             *
+             * \param[in] rec     Equilibration data of current region.
+             * \param[in] density Density calculator of current region.
+             * \param[in] rs      Calculator of dissolved gas-oil ratio.
+             * \param[in] rv      Calculator of vapourised oil-gas ratio.
+             * \param[in] pu      Summary of current active phases.
+             */
+            EquilReg(const EquilRecord& rec,
+                     const DensCalc&    density,
+                     std::shared_ptr<Miscibility::RsFunction> rs,
+                     std::shared_ptr<Miscibility::RsFunction> rv,
+                     const PhaseUsage&  pu)
+                : rec_    (rec)
+                , density_(density)
+                , rs_     (rs)
+                , rv_     (rv)
+                , pu_     (pu)
+            {
+            }
+
+            /**
+             * Type of density calculator.
+             */
+            typedef DensCalc CalcDensity;
+
+            /**
+             * Type of dissolved gas-oil ratio calculator.
+             */
+            typedef Miscibility::RsFunction CalcDissolution;
+
+            /**
+             * Type of vapourised oil-gas ratio calculator.
+             */
+            typedef Miscibility::RsFunction CalcEvaporation;
+
+            /**
+             * Datum depth in current region
+             */
+            double datum()    const { return this->rec_.main.depth; }
+
+            /**
+             * Pressure at datum depth in current region.
+             */
+            double pressure() const { return this->rec_.main.press; }
+
+            /**
+             * Depth of water-oil contact.
+             */
+            double zwoc()     const { return this->rec_.woc .depth; }
+
+            /**
+             * water-oil capillary pressure at water-oil contact.
+             *
+             * \return P_o - P_w at WOC.
+             */
+            double pcow_woc() const { return this->rec_.woc .press; }
+
+            /**
+             * Depth of gas-oil contact.
+             */
+            double zgoc()     const { return this->rec_.goc .depth; }
+
+            /**
+             * Gas-oil capillary pressure at gas-oil contact.
+             *
+             * \return P_g - P_o at GOC.
+             */
+            double pcgo_goc() const { return this->rec_.goc .press; }
+
+            /**
+             * Retrieve phase density calculator of current region.
+             */
+            const CalcDensity&
+            densityCalculator() const { return this->density_; }
+
+            /**
+             * Retrieve dissolved gas-oil ratio calculator of current
+             * region.
+             */
+            const CalcDissolution&
+            dissolutionCalculator() const { return *this->rs_; }
+
+            /**
+             * Retrieve vapourised oil-gas ratio calculator of current
+             * region.
+             */
+            const CalcEvaporation&
+            evaporationCalculator() const { return *this->rv_; }
+
+            /**
+             * Retrieve active fluid phase summary.
+             */
+            const PhaseUsage&
+            phaseUsage() const { return this->pu_; }
+
+        private:
+            EquilRecord rec_;     /**< Equilibration data */
+            DensCalc    density_; /**< Density calculator */
+            std::shared_ptr<Miscibility::RsFunction> rs_;      /**< RS calculator */
+            std::shared_ptr<Miscibility::RsFunction> rv_;      /**< RV calculator */
+            PhaseUsage  pu_;      /**< Active phase summary */
+        };
+
+
+
+        /// Functor for inverting capillary pressure function.
+        /// Function represented is
+        ///   f(s) = pc(s) - target_pc
+        struct PcEq
+        {
+            PcEq(const BlackoilPropertiesInterface& props,
+                 const int phase,
+                 const int cell,
+                 const double target_pc)
+                : props_(props),
+                  phase_(phase),
+                  cell_(cell),
+                  target_pc_(target_pc)
+            {
+                std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
+                std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
+            }
+            double operator()(double s) const
+            {
+                s_[phase_] = s;
+                props_.capPress(1, s_, &cell_, pc_, 0);
+                return pc_[phase_] - target_pc_;
+            }
+        private:
+            const BlackoilPropertiesInterface& props_;
+            const int phase_;
+            const int cell_;
+            const int target_pc_;
+            mutable double s_[BlackoilPhases::MaxNumPhases];
+            mutable double pc_[BlackoilPhases::MaxNumPhases];
+        };
+
+
+
+        /// Compute saturation of some phase corresponding to a given
+        /// capillary pressure.
+        inline double satFromPc(const BlackoilPropertiesInterface& props,
+                                const int phase,
+                                const int cell,
+                                const double target_pc,
+                                const bool increasing = false)
+        {
+            // Find minimum and maximum saturations.
+            double sminarr[BlackoilPhases::MaxNumPhases];
+            double smaxarr[BlackoilPhases::MaxNumPhases];
+            props.satRange(1, &cell, sminarr, smaxarr);
+            const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
+            const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
+
+            // Create the equation f(s) = pc(s) - target_pc
+            const PcEq f(props, phase, cell, target_pc);
+            const double f0 = f(s0);
+            const double f1 = f(s1);
+            if (f0 <= 0.0) {
+                return s0;
+            } else if (f1 > 0.0) {
+                return s1;
+            } else {
+                const int max_iter = 30;
+                const double tol = 1e-6;
+                int iter_used = -1;
+                typedef RegulaFalsi<ThrowOnError> ScalarSolver;
+                const double sol = ScalarSolver::solve(f, std::min(s0, s1), std::max(s0, s1), max_iter, tol, iter_used);
+                return sol;
+            }
+        }
+
+
+        /// Functor for inverting a sum of capillary pressure functions.
+        /// Function represented is
+        ///   f(s) = pc1(s) + pc2(1 - s) - target_pc
+        struct PcEqSum
+        {
+            PcEqSum(const BlackoilPropertiesInterface& props,
+                    const int phase1,
+                    const int phase2,
+                    const int cell,
+                    const double target_pc)
+                : props_(props),
+                  phase1_(phase1),
+                  phase2_(phase2),
+                  cell_(cell),
+                  target_pc_(target_pc)
+            {
+                std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
+                std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
+            }
+            double operator()(double s) const
+            {
+                s_[phase1_] = s;
+                s_[phase2_] = 1.0 - s;
+                props_.capPress(1, s_, &cell_, pc_, 0);
+                return pc_[phase1_] + pc_[phase2_] - target_pc_;
+            }
+        private:
+            const BlackoilPropertiesInterface& props_;
+            const int phase1_;
+            const int phase2_;
+            const int cell_;
+            const int target_pc_;
+            mutable double s_[BlackoilPhases::MaxNumPhases];
+            mutable double pc_[BlackoilPhases::MaxNumPhases];
+        };
+
+
+
+
+        /// Compute saturation of some phase corresponding to a given
+        /// capillary pressure, where the capillary pressure function
+        /// is given as a sum of two other functions.
+        inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
+                                      const int phase1,
+                                      const int phase2,
+                                      const int cell,
+                                      const double target_pc)
+        {
+            // Find minimum and maximum saturations.
+            double sminarr[BlackoilPhases::MaxNumPhases];
+            double smaxarr[BlackoilPhases::MaxNumPhases];
+            props.satRange(1, &cell, sminarr, smaxarr);
+            const double smin = sminarr[phase1];
+            const double smax = smaxarr[phase1];
+
+            // Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
+            const PcEqSum f(props, phase1, phase2, cell, target_pc);
+            const double f0 = f(smin);
+            const double f1 = f(smax);
+            if (f0 <= 0.0) {
+                return smin;
+            } else if (f1 > 0.0) {
+                return smax;
+            } else {
+                const int max_iter = 30;
+                const double tol = 1e-6;
+                int iter_used = -1;
+                typedef RegulaFalsi<ThrowOnError> ScalarSolver;
+                const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used);
+                return sol;
+            }
+        }
+
+    } // namespace Equil
+} // namespace Opm
+
+
+#endif // OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED

opm/core/simulator/initStateEquil.hpp

@@ -0,0 +1,673 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
+#define OPM_INITSTATEEQUIL_HEADER_INCLUDED
+
+#include <opm/core/simulator/EquilibrationHelpers.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/io/eclipse/EclipseGridParser.hpp>
+#include <opm/core/props/BlackoilPropertiesInterface.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/core/utility/RegionMapping.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/parser/eclipse/Utility/EquilWrapper.hpp>
+#include <opm/parser/eclipse/Utility/SimpleTable.hpp>
+
+#include <array>
+#include <cassert>
+#include <utility>
+#include <vector>
+
+/**
+ * \file
+ * Facilities for an ECLIPSE-style equilibration-based
+ * initialisation scheme (keyword 'EQUIL').
+ */
+struct UnstructuredGrid;
+
+namespace Opm
+{
+
+    /**
+     * Compute initial state by an equilibration procedure.
+     *
+     * The following state fields are modified:
+     *   pressure(),
+     *   saturation(),
+     *   surfacevol(),
+     *   gasoilratio(),
+     *   rv().
+     *
+     * \param[in] grid     Grid.
+     * \param[in] props    Property object, pvt and capillary properties are used.
+     * \param[in] deck     Simulation deck, used to obtain EQUIL and related data.
+     * \param[in] gravity  Acceleration of gravity, assumed to be in Z direction.
+     */
+    void initStateEquil(const UnstructuredGrid& grid,
+                        const BlackoilPropertiesInterface& props,
+                        const EclipseGridParser& deck,
+                        const double gravity,
+                        BlackoilState& state);
+                        
+    
+    void initStateEquil(const UnstructuredGrid& grid,
+                        const BlackoilPropertiesInterface& props,
+                        const Opm::DeckConstPtr newParserDeck,
+                        const double gravity,
+                        BlackoilState& state);
+
+
+    /**
+     * Types and routines that collectively implement a basic
+     * ECLIPSE-style equilibration-based initialisation scheme.
+     *
+     * This namespace is intentionally nested to avoid name clashes
+     * with other parts of OPM.
+     */
+    namespace Equil {
+
+        /**
+         * Compute initial phase pressures by means of equilibration.
+         *
+         * This function uses the information contained in an
+         * equilibration record (i.e., depths and pressurs) as well as
+         * a density calculator and related data to vertically
+         * integrate the phase pressure ODE
+         * \f[
+         * \frac{\mathrm{d}p_{\alpha}}{\mathrm{d}z} =
+         * \rho_{\alpha}(z,p_{\alpha})\cdot g
+         * \f]
+         * in which \f$\rho_{\alpha}$ denotes the fluid density of
+         * fluid phase \f$\alpha\f$, \f$p_{\alpha}\f$ is the
+         * corresponding phase pressure, \f$z\f$ is the depth and
+         * \f$g\f$ is the acceleration due to gravity (assumed
+         * directed downwords, in the positive \f$z\f$ direction).
+         *
+         * \tparam Region Type of an equilibration region information
+         *                base.  Typically an instance of the EquilReg
+         *                class template.
+         *
+         * \tparam CellRange Type of cell range that demarcates the
+         *                cells pertaining to the current
+         *                equilibration region.  Must implement
+         *                methods begin() and end() to bound the range
+         *                as well as provide an inner type,
+         *                const_iterator, to traverse the range.
+         *
+         * \param[in] G     Grid.
+         * \param[in] reg   Current equilibration region.
+         * \param[in] cells Range that spans the cells of the current
+         *                  equilibration region.
+         * \param[in] grav  Acceleration of gravity.
+         *
+         * \return Phase pressures, one vector for each active phase,
+         * of pressure values in each cell in the current
+         * equilibration region.
+         */
+        template <class Region, class CellRange>
+        std::vector< std::vector<double> >
+        phasePressures(const UnstructuredGrid& G,
+                       const Region&           reg,
+                       const CellRange&        cells,
+                       const double            grav = unit::gravity);
+
+
+
+        /**
+         * Compute initial phase saturations by means of equilibration.
+         *
+         * \tparam Region Type of an equilibration region information
+         *                base.  Typically an instance of the EquilReg
+         *                class template.
+         *
+         * \tparam CellRange Type of cell range that demarcates the
+         *                cells pertaining to the current
+         *                equilibration region.  Must implement
+         *                methods begin() and end() to bound the range
+         *                as well as provide an inner type,
+         *                const_iterator, to traverse the range.
+         *
+         * \param[in] reg             Current equilibration region.
+         * \param[in] cells           Range that spans the cells of the current
+         *                            equilibration region.
+         * \param[in] props           Property object, needed for capillary functions.
+         * \param[in] phase_pressures Phase pressures, one vector for each active phase,
+         *                            of pressure values in each cell in the current
+         *                            equilibration region.
+         * \return                    Phase saturations, one vector for each phase, each containing
+         *                            one saturation value per cell in the region.
+         */
+        template <class Region, class CellRange>
+        std::vector< std::vector<double> >
+        phaseSaturations(const Region&           reg,
+                         const CellRange&        cells,
+                         const BlackoilPropertiesInterface& props,
+                         const std::vector< std::vector<double> >& phase_pressures);
+
+
+
+        /**
+         * Compute initial Rs values.
+         *
+         * \tparam CellRangeType Type of cell range that demarcates the
+         *                cells pertaining to the current
+         *                equilibration region.  Must implement
+         *                methods begin() and end() to bound the range
+         *                as well as provide an inner type,
+         *                const_iterator, to traverse the range.
+         *
+         * \param[in] grid            Grid.
+         * \param[in] cells           Range that spans the cells of the current
+         *                            equilibration region.
+         * \param[in] oil_pressure    Oil pressure for each cell in range.
+         * \param[in] rs_func         Rs as function of pressure and depth.
+         * \return                    Rs values, one for each cell in the 'cells' range.
+         */
+        template <class CellRangeType>
+        std::vector<double> computeRs(const UnstructuredGrid& grid,
+                                      const CellRangeType& cells,
+                                      const std::vector<double> oil_pressure,
+                                      const Miscibility::RsFunction& rs_func,
+                                      const std::vector<double> gas_saturation);
+
+        namespace DeckDependent {
+
+            inline
+            std::vector<EquilRecord>
+            getEquil(const EclipseGridParser& deck)
+            {
+                if (deck.hasField("EQUIL")) {
+                    const EQUIL& eql = deck.getEQUIL();
+
+                    typedef std::vector<EquilLine>::size_type sz_t;
+                    const sz_t nrec = eql.equil.size();
+
+                    std::vector<EquilRecord> ret;
+                    ret.reserve(nrec);
+                    for (sz_t r = 0; r < nrec; ++r) {
+                        const EquilLine& rec = eql.equil[r];
+
+                        EquilRecord record =
+                            {
+                                { rec.datum_depth_             ,
+                                  rec.datum_depth_pressure_    }
+                                ,
+                                { rec.water_oil_contact_depth_ ,
+                                  rec.oil_water_cap_pressure_  }
+                                ,
+                                { rec.gas_oil_contact_depth_   ,
+                                  rec.gas_oil_cap_pressure_    }
+                                ,
+                                rec.live_oil_table_index_
+                                ,
+                                rec.wet_gas_table_index_
+                                ,
+                                rec.N_
+                            };
+                        if (record.N != 0) {
+                            OPM_THROW(std::domain_error,
+                              "kw EQUIL, item 9: Only N=0 supported.");
+                        }
+                        ret.push_back(record);
+                    }
+
+                    return ret;
+                }
+                else {
+                    OPM_THROW(std::domain_error,
+                              "Deck does not provide equilibration data.");
+                }
+            }
+
+            inline
+            std::vector<EquilRecord>
+            getEquil(const Opm::DeckConstPtr newParserDeck)
+            {
+                if (newParserDeck->hasKeyword("EQUIL")) {
+                
+                    Opm::EquilWrapper eql(newParserDeck->getKeyword("EQUIL"));
+
+                    const int nrec = eql.numRegions();
+
+                    std::vector<EquilRecord> ret;
+                    ret.reserve(nrec);
+                    for (int r = 0; r < nrec; ++r) {
+
+                        EquilRecord record =
+                            {
+                                { eql.datumDepth(r)                        ,
+                                  eql.datumDepthPressure(r)                }
+                                ,
+                                { eql.waterOilContactDepth(r)              ,
+                                  eql.waterOilContactCapillaryPressure(r)  }
+                                ,
+                                { eql.gasOilContactDepth(r)                ,
+                                  eql.gasOilContactCapillaryPressure(r)    }
+                                ,
+                                eql.liveOilInitProceedure(r)
+                                ,
+                                eql.wetGasInitProceedure(r)
+                                ,
+                                eql.initializationTargetAccuracy(r)
+                            };
+                        if (record.N != 0) {
+                            OPM_THROW(std::domain_error,
+                              "kw EQUIL, item 9: Only N=0 supported.");
+                        }
+                        ret.push_back(record);
+                    }
+
+                    return ret;
+                }
+                else {
+                    OPM_THROW(std::domain_error,
+                              "Deck does not provide equilibration data.");
+                }
+            }
+
+
+            inline
+            std::vector<int>
+            equilnum(const EclipseGridParser& deck,
+                     const UnstructuredGrid&  G   )
+            {
+                std::vector<int> eqlnum;
+                if (deck.hasField("EQLNUM")) {        
+                    const std::vector<int>& e = deck.getIntegerValue("EQLNUM");
+                    eqlnum.reserve(e.size());
+                    std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
+                                   std::bind2nd(std::minus<int>(), 1));
+                }
+                else {
+                    // No explicit equilibration region.
+                    // All cells in region zero.
+                    eqlnum.assign(G.number_of_cells, 0);
+                }
+
+                return eqlnum;
+            }
+
+
+            inline
+            std::vector<int>
+            equilnum(const Opm::DeckConstPtr newParserDeck,
+                     const UnstructuredGrid&  G   )
+            {
+                std::vector<int> eqlnum;
+                if (newParserDeck->hasKeyword("EQLNUM")) {        
+                    const std::vector<int>& e = newParserDeck->getKeyword("EQLNUM")->getIntData();
+                    eqlnum.reserve(e.size());
+                    std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
+                                   std::bind2nd(std::minus<int>(), 1));
+                }
+                else {
+                    // No explicit equilibration region.
+                    // All cells in region zero.
+                    eqlnum.assign(G.number_of_cells, 0);
+                }
+
+                return eqlnum;
+            }
+
+
+            template <class InputDeck>
+            class InitialStateComputer;
+
+            template <>
+            class InitialStateComputer<Opm::EclipseGridParser> {
+            public:
+                InitialStateComputer(const BlackoilPropertiesInterface& props,
+                                     const EclipseGridParser&           deck ,
+                                     const UnstructuredGrid&            G    ,
+                                     const double                       grav = unit::gravity)
+                    : pp_(props.numPhases(),
+                          std::vector<double>(G.number_of_cells)),
+                      sat_(props.numPhases(),
+                          std::vector<double>(G.number_of_cells)),
+                      rs_(G.number_of_cells),
+                      rv_(G.number_of_cells)
+                {
+                    // Get the equilibration records.
+                    const std::vector<EquilRecord> rec = getEquil(deck);
+
+                    // Create (inverse) region mapping.
+                    const RegionMapping<> eqlmap(equilnum(deck, G));
+
+                    // Create Rs functions.
+                    rs_func_.reserve(rec.size());
+                    if (deck.hasField("DISGAS")) {                    
+                        for (size_t i = 0; i < rec.size(); ++i) {
+                            const int cell = *(eqlmap.cells(i).begin());                   
+                            if (rec[i].live_oil_table_index > 0) {
+                                if (deck.hasField("RSVD")) { 
+                                    // TODO When this kw is actually parsed, also check for proper number of available tables
+                                    // For now, just use dummy ...
+                                    std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
+                                    std::vector<double> rs; rs.push_back(0.0); rs.push_back(100.0);
+                                    rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, depth, rs));
+                                } else {
+                                    OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
+                                }
+                            } else {
+                                if (rec[i].goc.depth != rec[i].main.depth) {
+                                    OPM_THROW(std::runtime_error,
+                                              "Cannot initialise: when no explicit RSVD table is given, \n"
+                                              "datum depth must be at the gas-oil-contact. "
+                                              "In EQUIL region " << (i + 1) << "  (counting from 1), this does not hold.");
+                                }
+                                const double p_contact = rec[i].main.press;
+                                rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact));
+                            }
+                        }
+                    } else {
+                        for (size_t i = 0; i < rec.size(); ++i) {
+                            rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
+                        }
+                    }                    
+
+                    rv_func_.reserve(rec.size());
+                    if (deck.hasField("VAPOIL")) {                    
+                        for (size_t i = 0; i < rec.size(); ++i) {
+                            const int cell = *(eqlmap.cells(i).begin());                   
+                            if (rec[i].wet_gas_table_index > 0) {
+                                if (deck.hasField("RVVD")) { 
+                                    // TODO When this kw is actually parsed, also check for proper number of available tables
+                                    // For now, just use dummy ...
+                                    std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
+                                    std::vector<double> rv; rv.push_back(0.0); rv.push_back(0.0001);
+                                    rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, depth, rv));
+                                } else {
+                                    OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
+                                }
+                            } else {
+                                if (rec[i].goc.depth != rec[i].main.depth) {
+                                    OPM_THROW(std::runtime_error,
+                                              "Cannot initialise: when no explicit RVVD table is given, \n"
+                                              "datum depth must be at the gas-oil-contact. "
+                                              "In EQUIL region " << (i + 1) << "  (counting from 1), this does not hold.");
+                                }
+                                const double p_contact = rec[i].main.press + rec[i].goc.press;
+                                rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact));
+                            }
+                        }
+                    } else {
+                        for (size_t i = 0; i < rec.size(); ++i) {
+                            rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
+                        }
+                    }
+                    
+                    // Compute pressures, saturations, rs and rv factors.
+                    calcPressSatRsRv(eqlmap, rec, props, G, grav);
+
+                    // Modify oil pressure in no-oil regions so that the pressures of present phases can
+                    // be recovered from the oil pressure and capillary relations.
+                }
+
+                typedef std::vector<double> Vec;
+                typedef std::vector<Vec>    PVec; // One per phase.
+
+                const PVec& press() const { return pp_; }
+                const PVec& saturation() const { return sat_; }
+                const Vec& rs() const { return rs_; }
+                const Vec& rv() const { return rv_; }
+
+            private:
+                typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
+                typedef EquilReg<RhoCalc> EqReg;
+
+                std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
+                std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
+
+                PVec pp_;
+                PVec sat_;
+                Vec rs_;
+                Vec rv_;
+
+                template <class RMap>
+                void
+                calcPressSatRsRv(const RMap&                             reg  ,
+                                 const std::vector< EquilRecord >&       rec  ,
+                                 const Opm::BlackoilPropertiesInterface& props,
+                                 const UnstructuredGrid&                 G    ,
+                                 const double grav)
+                {
+                    typedef Miscibility::NoMixing NoMix;
+
+                    for (typename RMap::RegionId
+                             r = 0, nr = reg.numRegions();
+                         r < nr; ++r)
+                    {
+                        const typename RMap::CellRange cells = reg.cells(r);
+
+                        const int repcell = *cells.begin();
+                        const RhoCalc calc(props, repcell);
+                        const EqReg eqreg(rec[r], calc,
+                                          rs_func_[r], rv_func_[r],
+                                          props.phaseUsage());
+                                          
+                        PVec press = phasePressures(G, eqreg, cells, grav);
+                        
+                        const PVec sat = phaseSaturations(eqreg, cells, props, press);
+                        
+                        const int np = props.numPhases();
+                        for (int p = 0; p < np; ++p) {
+                            copyFromRegion(press[p], cells, pp_[p]);
+                            copyFromRegion(sat[p], cells, sat_[p]);
+                        }
+                        if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
+                            && props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
+                            const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+                            const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+                            const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
+                            const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
+                            copyFromRegion(rs, cells, rs_);
+                            copyFromRegion(rv, cells, rv_);
+                        }
+                    }
+                }
+
+                template <class CellRangeType>
+                void copyFromRegion(const Vec& source,
+                                    const CellRangeType& cells,
+                                    Vec& destination)
+                {
+                    auto s = source.begin();
+                    auto c = cells.begin();
+                    const auto e = cells.end();
+                    for (; c != e; ++c, ++s) {
+                        destination[*c] = *s;
+                    }
+                }
+
+            };
+
+
+            template <>
+            class InitialStateComputer<Opm::DeckConstPtr> {
+            public:
+                InitialStateComputer(const BlackoilPropertiesInterface& props,
+                                     const Opm::DeckConstPtr            newParserDeck,
+                                     const UnstructuredGrid&            G    ,
+                                     const double                       grav = unit::gravity)
+                    : pp_(props.numPhases(),
+                          std::vector<double>(G.number_of_cells)),
+                      sat_(props.numPhases(),
+                          std::vector<double>(G.number_of_cells)),
+                      rs_(G.number_of_cells),
+                      rv_(G.number_of_cells)
+                {
+                    // Get the equilibration records.
+                    const std::vector<EquilRecord> rec = getEquil(newParserDeck);
+
+                    // Create (inverse) region mapping.
+                    const RegionMapping<> eqlmap(equilnum(newParserDeck, G)); 
+
+                    // Create Rs functions.
+                    rs_func_.reserve(rec.size());
+                    if (newParserDeck->hasKeyword("DISGAS")) {                    
+                        for (size_t i = 0; i < rec.size(); ++i) {
+                            const int cell = *(eqlmap.cells(i).begin());                   
+                            if (rec[i].live_oil_table_index > 0) {
+                                if (newParserDeck->hasKeyword("RSVD") && rec[i].live_oil_table_index <= newParserDeck->getKeyword("RSVD")->size()) { 
+                                    Opm::SimpleTable rsvd(newParserDeck->getKeyword("RSVD"),std::vector<std::string>{"vd", "rs"},rec[i].live_oil_table_index-1);                                
+                                    std::vector<double> vd(rsvd.getColumn("vd"));
+                                    std::vector<double> rs(rsvd.getColumn("rs"));
+                                    rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, vd, rs));
+                                } else {
+                                    OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
+                                }
+                            } else {
+                                if (rec[i].goc.depth != rec[i].main.depth) {
+                                    OPM_THROW(std::runtime_error,
+                                              "Cannot initialise: when no explicit RSVD table is given, \n"
+                                              "datum depth must be at the gas-oil-contact. "
+                                              "In EQUIL region " << (i + 1) << "  (counting from 1), this does not hold.");
+                                }
+                                const double p_contact = rec[i].main.press;
+                                rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact));
+                            }
+                        }
+                    } else {
+                        for (size_t i = 0; i < rec.size(); ++i) {
+                            rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
+                        }
+                    }                    
+
+                    rv_func_.reserve(rec.size());
+                    if (newParserDeck->hasKeyword("VAPOIL")) {                    
+                        for (size_t i = 0; i < rec.size(); ++i) {
+                            const int cell = *(eqlmap.cells(i).begin());                   
+                            if (rec[i].wet_gas_table_index > 0) {
+                                if (newParserDeck->hasKeyword("RVVD") && rec[i].wet_gas_table_index <= newParserDeck->getKeyword("RVVD")->size()) { 
+                                    Opm::SimpleTable rvvd(newParserDeck->getKeyword("RVVD"),std::vector<std::string>{"vd", "rv"},rec[i].wet_gas_table_index-1);                                
+                                    std::vector<double> vd(rvvd.getColumn("vd"));
+                                    std::vector<double> rv(rvvd.getColumn("rv"));
+                                    rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, vd, rv));
+                                } else {
+                                    OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
+                                }
+                            } else {
+                                if (rec[i].goc.depth != rec[i].main.depth) {
+                                    OPM_THROW(std::runtime_error,
+                                              "Cannot initialise: when no explicit RVVD table is given, \n"
+                                              "datum depth must be at the gas-oil-contact. "
+                                              "In EQUIL region " << (i + 1) << "  (counting from 1), this does not hold.");
+                                }
+                                const double p_contact = rec[i].main.press + rec[i].goc.press;
+                                rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact));
+                            }
+                        }
+                    } else {
+                        for (size_t i = 0; i < rec.size(); ++i) {
+                            rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
+                        }
+                    }
+
+                    // Compute pressures, saturations, rs and rv factors.
+                    calcPressSatRsRv(eqlmap, rec, props, G, grav);
+
+                    // Modify oil pressure in no-oil regions so that the pressures of present phases can
+                    // be recovered from the oil pressure and capillary relations.
+                }
+
+                typedef std::vector<double> Vec;
+                typedef std::vector<Vec>    PVec; // One per phase.
+
+                const PVec& press() const { return pp_; }
+                const PVec& saturation() const { return sat_; }
+                const Vec& rs() const { return rs_; }
+                const Vec& rv() const { return rv_; }
+
+            private:
+                typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
+                typedef EquilReg<RhoCalc> EqReg;
+
+                std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
+                std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
+
+                PVec pp_;
+                PVec sat_;
+                Vec rs_;
+                Vec rv_;
+
+                template <class RMap>
+                void
+                calcPressSatRsRv(const RMap&                             reg  ,
+                                 const std::vector< EquilRecord >&       rec  ,
+                                 const Opm::BlackoilPropertiesInterface& props,
+                                 const UnstructuredGrid&                 G    ,
+                                 const double grav)
+                {
+                    typedef Miscibility::NoMixing NoMix;
+
+                    for (typename RMap::RegionId
+                             r = 0, nr = reg.numRegions();
+                         r < nr; ++r)
+                    {
+                        const typename RMap::CellRange cells = reg.cells(r);
+
+                        const int repcell = *cells.begin();
+                        const RhoCalc calc(props, repcell);
+                        const EqReg eqreg(rec[r], calc,
+                                          rs_func_[r], rv_func_[r],
+                                          props.phaseUsage());
+                   
+                        PVec press = phasePressures(G, eqreg, cells, grav);
+
+                        const PVec sat = phaseSaturations(eqreg, cells, props, press);
+
+                        const int np = props.numPhases();
+                        for (int p = 0; p < np; ++p) {
+                            copyFromRegion(press[p], cells, pp_[p]);
+                            copyFromRegion(sat[p], cells, sat_[p]);
+                        }
+                        if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
+                            && props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
+                            const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+                            const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+                            const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
+                            const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
+                            copyFromRegion(rs, cells, rs_);
+                            copyFromRegion(rv, cells, rv_);
+                        }
+                    }
+                }
+
+                template <class CellRangeType>
+                void copyFromRegion(const Vec& source,
+                                    const CellRangeType& cells,
+                                    Vec& destination)
+                {
+                    auto s = source.begin();
+                    auto c = cells.begin();
+                    const auto e = cells.end();
+                    for (; c != e; ++c, ++s) {
+                        destination[*c] = *s;
+                    }
+                }
+
+            };
+        } // namespace DeckDependent
+    } // namespace Equil
+} // namespace Opm
+
+#include <opm/core/simulator/initStateEquil_impl.hpp>
+
+#endif // OPM_INITSTATEEQUIL_HEADER_INCLUDED

opm/core/simulator/initStateEquil_impl.hpp

@@ -0,0 +1,781 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
+#define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
+
+#include <opm/core/grid.h>
+#include <opm/core/props/BlackoilPhases.hpp>
+
+#include <cassert>
+#include <cmath>
+#include <functional>
+#include <vector>
+
+namespace Opm
+{
+    namespace Details {
+        template <class RHS>
+        class RK4IVP {
+        public:
+            RK4IVP(const RHS&                  f   ,
+                   const std::array<double,2>& span,
+                   const double                y0  ,
+                   const int                   N   )
+                : N_(N)
+                , span_(span)
+            {
+                const double h  = stepsize();
+                const double h2 = h / 2;
+                const double h6 = h / 6;
+
+                y_.reserve(N + 1);
+                f_.reserve(N + 1);
+
+                y_.push_back(y0);
+                f_.push_back(f(span_[0], y0));
+
+                for (int i = 0; i < N; ++i) {
+                    const double x  = span_[0] + i*h;
+                    const double y  = y_.back();
+
+                    const double k1 = f_[i];
+                    const double k2 = f(x + h2, y + h2*k1);
+                    const double k3 = f(x + h2, y + h2*k2);
+                    const double k4 = f(x + h , y + h*k3);
+
+                    y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4));
+                    f_.push_back(f(x + h, y_.back()));
+                }
+
+                assert (y_.size() == std::vector<double>::size_type(N + 1));
+            }
+
+            double
+            operator()(const double x) const
+            {
+                // Dense output (O(h**3)) according to Shampine
+                // (Hermite interpolation)
+                const double h = stepsize();
+                int          i = (x - span_[0]) / h;
+                const double t = (x - (span_[0] + i*h)) / h;
+
+                // Crude handling of evaluation point outside "span_";
+                if (i  <  0) { i = 0;      }
+                if (N_ <= i) { i = N_ - 1; }
+
+                const double y0 = y_[i], y1 = y_[i + 1];
+                const double f0 = f_[i], f1 = f_[i + 1];
+
+                double u = (1 - 2*t) * (y1 - y0);
+                u += h * ((t - 1)*f0 + t*f1);
+                u *= t * (t - 1);
+                u += (1 - t)*y0 + t*y1;
+
+                return u;
+            }
+
+        private:
+            int                  N_;
+            std::array<double,2> span_;
+            std::vector<double>  y_;
+            std::vector<double>  f_;
+
+            double
+            stepsize() const { return (span_[1] - span_[0]) / N_; }
+        };
+
+        namespace PhasePressODE {
+            template <class Density>
+            class Water {
+            public:
+                Water(const Density& rho,
+                      const int      np,
+                      const int      ix,
+                      const double   norm_grav)
+                    : rho_(rho)
+                    , svol_(np, 0)
+                    , ix_(ix)
+                    , g_(norm_grav)
+                {
+                    svol_[ix_] = 1.0;
+                }
+
+                double
+                operator()(const double /* depth */,
+                           const double press) const
+                {
+                    return this->density(press) * g_;
+                }
+
+            private:
+                const Density&      rho_;
+                std::vector<double> svol_;
+                const int           ix_;
+                const double        g_;
+
+                double
+                density(const double press) const
+                {
+                    const std::vector<double>& rho = rho_(press, svol_);
+
+                    return rho[ix_];
+                }
+            };
+
+            template <class Density, class RS>
+            class Oil {
+            public:
+                Oil(const Density& rho,
+                    const RS&      rs,
+                    const int      np,
+                    const int      oix,
+                    const int      gix,
+                    const double   norm_grav)
+                    : rho_(rho)
+                    , rs_(rs)
+                    , svol_(np, 0)
+                    , oix_(oix)
+                    , gix_(gix)
+                    , g_(norm_grav)
+                {
+                    svol_[oix_] = 1.0;
+                }
+
+                double
+                operator()(const double depth,
+                           const double press) const
+                {
+                    return this->density(depth, press) * g_;
+                }
+
+            private:
+                const Density&              rho_;
+                const RS&                   rs_;
+                mutable std::vector<double> svol_;
+                const int                   oix_;
+                const int                   gix_;
+                const double                g_;
+
+                double
+                density(const double depth,
+                        const double press) const
+                {
+                    if (gix_ >= 0) {
+                        svol_[gix_] = rs_(depth, press);
+                    }
+
+                    const std::vector<double>& rho = rho_(press, svol_);
+                    return rho[oix_];
+                }
+            };
+
+            template <class Density, class RV>
+            class Gas {
+            public:
+                Gas(const Density& rho,
+                    const RV&      rv,
+                    const int      np,
+                    const int      gix,
+                    const int      oix,
+                    const double   norm_grav)
+                    : rho_(rho)
+                    , rv_(rv)
+                    , svol_(np, 0)
+                    , gix_(gix)
+                    , oix_(oix)
+                    , g_(norm_grav)
+                {
+                    svol_[gix_] = 1.0;
+                }
+
+                double
+                operator()(const double depth,
+                           const double press) const
+                {
+                    return this->density(depth, press) * g_;
+                }
+
+            private:
+                const Density&              rho_;
+                const RV&                   rv_;
+                mutable std::vector<double> svol_;
+                const int                   gix_;
+                const int                   oix_;
+                const double                g_;
+
+                double
+                density(const double depth,
+                        const double press) const
+                {
+                    if (oix_ >= 0) {
+                        svol_[oix_] = rv_(depth, press);
+                    }
+
+                    const std::vector<double>& rho = rho_(press, svol_);
+                    return rho[gix_];
+                }
+            };
+        } // namespace PhasePressODE
+
+        namespace PhaseUsed {
+            inline bool
+            water(const PhaseUsage& pu)
+            {
+                return bool(pu.phase_used[ Opm::BlackoilPhases::Aqua ]);
+            }
+
+            inline bool
+            oil(const PhaseUsage& pu)
+            {
+                return bool(pu.phase_used[ Opm::BlackoilPhases::Liquid ]);
+            }
+
+            inline bool
+            gas(const PhaseUsage& pu)
+            {
+                return bool(pu.phase_used[ Opm::BlackoilPhases::Vapour ]);
+            }
+        } // namespace PhaseUsed
+
+        namespace PhaseIndex {
+            inline int
+            water(const PhaseUsage& pu)
+            {
+                int i = -1;
+                if (PhaseUsed::water(pu)) {
+                    i = pu.phase_pos[ Opm::BlackoilPhases::Aqua ];
+                }
+
+                return i;
+            }
+
+            inline int
+            oil(const PhaseUsage& pu)
+            {
+                int i = -1;
+                if (PhaseUsed::oil(pu)) {
+                    i = pu.phase_pos[ Opm::BlackoilPhases::Liquid ];
+                }
+
+                return i;
+            }
+
+            inline int
+            gas(const PhaseUsage& pu)
+            {
+                int i = -1;
+                if (PhaseUsed::gas(pu)) {
+                    i = pu.phase_pos[ Opm::BlackoilPhases::Vapour ];
+                }
+
+                return i;
+            }
+        } // namespace PhaseIndex
+
+        namespace PhasePressure {
+            template <class PressFunction,
+                      class CellRange>
+            void
+            assign(const UnstructuredGrid&             G    ,
+                   const std::array<PressFunction, 2>& f    ,
+                   const double                        split,
+                   const CellRange&                    cells,
+                   std::vector<double>&                p    )
+            {
+                const int nd = G.dimensions;
+
+                enum { up = 0, down = 1 };
+
+                std::vector<double>::size_type c = 0;
+                for (typename CellRange::const_iterator
+                         ci = cells.begin(), ce = cells.end();
+                     ci != ce; ++ci, ++c)
+                {
+                    assert (c < p.size());
+
+                    const double z = G.cell_centroids[(*ci)*nd + (nd - 1)];
+                    p[c] = (z < split) ? f[up](z) : f[down](z);
+                }
+            }
+
+            template <class Region,
+                      class CellRange>
+            void
+            water(const UnstructuredGrid&     G     ,
+                  const Region&               reg   ,
+                  const std::array<double,2>& span  ,
+                  const double                grav  ,
+                  const double                po_woc,
+                  const CellRange&            cells ,
+                  std::vector<double>&        press )
+            {
+                using PhasePressODE::Water;
+                typedef Water<typename Region::CalcDensity> ODE;
+
+                const PhaseUsage& pu = reg.phaseUsage();
+
+                const int wix = PhaseIndex::water(pu);
+                ODE drho(reg.densityCalculator(), pu.num_phases, wix, grav);
+
+                const double z0 = reg.zwoc();
+                const double p0 = po_woc - reg.pcow_woc(); // Pcow = Po - Pw
+
+                std::array<double,2> up   = {{ z0, span[0] }};
+                std::array<double,2> down = {{ z0, span[1] }};
+
+                typedef Details::RK4IVP<ODE> WPress;
+                std::array<WPress,2> wpress = {
+                    {
+                        WPress(drho, up  , p0, 100)
+                        ,
+                        WPress(drho, down, p0, 100)
+                    }
+                };
+
+                assign(G, wpress, z0, cells, press);
+            }
+
+            template <class Region,
+                      class CellRange>
+            void
+            oil(const UnstructuredGrid&     G     ,
+                const Region&               reg   ,
+                const std::array<double,2>& span  ,
+                const double                grav  ,
+                const CellRange&            cells ,
+                std::vector<double>&        press ,
+                double&                     po_woc,
+                double&                     po_goc)
+            {
+                using PhasePressODE::Oil;
+                typedef Oil<typename Region::CalcDensity,
+                            typename Region::CalcDissolution> ODE;
+
+                const PhaseUsage& pu = reg.phaseUsage();
+
+                const int oix = PhaseIndex::oil(pu);
+                const int gix = PhaseIndex::gas(pu);
+                ODE drho(reg.densityCalculator(),
+                         reg.dissolutionCalculator(),
+                         pu.num_phases, oix, gix, grav);
+
+                const double z0 = reg.datum();
+                const double p0 = reg.pressure();
+
+                std::array<double,2> up   = {{ z0, span[0] }};
+                std::array<double,2> down = {{ z0, span[1] }};
+
+                typedef Details::RK4IVP<ODE> OPress;
+                std::array<OPress,2> opress = {
+                    {
+                        OPress(drho, up  , p0, 100)
+                        ,
+                        OPress(drho, down, p0, 100)
+                    }
+                };
+
+                assign(G, opress, z0, cells, press);
+
+                po_woc = -std::numeric_limits<double>::max();
+                po_goc = -std::numeric_limits<double>::max();
+
+                const double woc = reg.zwoc();
+                // Compute Oil pressure at WOC
+                if      (z0 > woc) { po_woc = opress[0](woc); } // WOC above datum
+                else if (z0 < woc) { po_woc = opress[1](woc); } // WOC below datum
+                else               { po_woc = p0;             } // WOC *at*  datum
+
+                const double goc = reg.zgoc();
+                // Compute Oil pressure at GOC
+                if      (z0 > goc) { po_goc = opress[0](goc); } // GOC above datum
+                else if (z0 < goc) { po_goc = opress[1](goc); } // GOC below datum
+                else               { po_goc = p0;             } // GOC *at*  datum
+            }
+
+            template <class Region,
+                      class CellRange>
+            void
+            gas(const UnstructuredGrid&     G     ,
+                const Region&               reg   ,
+                const std::array<double,2>& span  ,
+                const double                grav  ,
+                const double                po_goc,
+                const CellRange&            cells ,
+                std::vector<double>&        press )
+            {
+                using PhasePressODE::Gas;
+                typedef Gas<typename Region::CalcDensity,
+                            typename Region::CalcEvaporation> ODE;
+
+                const PhaseUsage& pu = reg.phaseUsage();
+
+                const int gix = PhaseIndex::gas(pu);
+                const int oix = PhaseIndex::oil(pu);
+
+                ODE drho(reg.densityCalculator(),
+                         reg.evaporationCalculator(),
+                         pu.num_phases, gix, oix, grav);
+
+                const double z0 = reg.zgoc();
+                const double p0 = po_goc + reg.pcgo_goc(); // Pcog = Pg - Po
+
+                std::array<double,2> up   = {{ z0, span[0] }};
+                std::array<double,2> down = {{ z0, span[1] }};
+
+                typedef Details::RK4IVP<ODE> GPress;
+                std::array<GPress,2> gpress = {
+                    {
+                        GPress(drho, up  , p0, 100)
+                        ,
+                        GPress(drho, down, p0, 100)
+                    }
+                };
+
+                assign(G, gpress, z0, cells, press);
+            }
+        } // namespace PhasePressure
+
+        template <class Region,
+                  class CellRange>
+        void
+        equilibrateOWG(const UnstructuredGrid&             G,
+                       const Region&                       reg,
+                       const double                        grav,
+                       const std::array<double,2>&         span,
+                       const CellRange&                    cells,
+                       std::vector< std::vector<double> >& press)
+        {
+            const PhaseUsage& pu = reg.phaseUsage();
+
+            double po_woc = -1, po_goc = -1;
+            if (PhaseUsed::oil(pu)) {
+                const int oix = PhaseIndex::oil(pu);
+                PhasePressure::oil(G, reg, span, grav, cells,
+                                   press[ oix ], po_woc, po_goc);
+            }
+
+            if (PhaseUsed::water(pu)) {
+                const int wix = PhaseIndex::water(pu);
+                PhasePressure::water(G, reg, span, grav, po_woc,
+                                     cells, press[ wix ]);
+            }
+
+            if (PhaseUsed::gas(pu)) {
+                const int gix = PhaseIndex::gas(pu);
+                PhasePressure::gas(G, reg, span, grav, po_goc,
+                                   cells, press[ gix ]);
+            }
+        }
+    } // namespace Details
+
+
+    namespace Equil {
+
+
+        template <class Region,
+                  class CellRange>
+        std::vector< std::vector<double> >
+        phasePressures(const UnstructuredGrid& G,
+                       const Region&           reg,
+                       const CellRange&        cells,
+                       const double            grav)
+        {
+            std::array<double,2> span =
+                {{  std::numeric_limits<double>::max() ,
+                   -std::numeric_limits<double>::max() }}; // Symm. about 0.
+
+            int ncell = 0;
+            {
+                // This code is only supported in three space dimensions
+                assert (G.dimensions == 3);
+
+                const int nd = G.dimensions;
+
+                // Define short-name aliases to reduce visual clutter.
+                const double* const nc  = & G.node_coordinates[0];
+
+                const int*    const cfp = & G.cell_facepos[0];
+                const int*    const cf  = & G.cell_faces[0];
+
+                const int*    const fnp = & G.face_nodepos[0];
+                const int*    const fn  = & G.face_nodes[0];
+
+                // Define vertical span as
+                //
+                //   [minimum(node depth(cells)), maximum(node depth(cells))]
+                //
+                // Note: We use a sledgehammer approach--looping all
+                // the nodes of all the faces of all the 'cells'--to
+                // compute those bounds.  This necessarily entails
+                // visiting some nodes (and faces) multiple times.
+                //
+                // Note: The implementation of 'RK4IVP<>' implicitly
+                // imposes the requirement that cell centroids are all
+                // within this vertical span.  That requirement is not
+                // checked.
+                for (typename CellRange::const_iterator
+                         ci = cells.begin(), ce = cells.end();
+                     ci != ce; ++ci, ++ncell)
+                {
+                    for (const int
+                             *fi = & cf[ cfp[*ci + 0] ],
+                             *fe = & cf[ cfp[*ci + 1] ];
+                         fi != fe; ++fi)
+                    {
+                        for (const int
+                                 *i = & fn[ fnp[*fi + 0] ],
+                                 *e = & fn[ fnp[*fi + 1] ];
+                             i != e; ++i)
+                        {
+                            const double z = nc[(*i)*nd + (nd - 1)];
+
+                            if (z < span[0]) { span[0] = z; }
+                            if (z > span[1]) { span[1] = z; }
+                        }
+                    }
+                }
+            }
+
+            const int np = reg.phaseUsage().num_phases;
+
+            typedef std::vector<double> pval;
+            std::vector<pval> press(np, pval(ncell, 0.0));
+
+            const double z0   = reg.datum();
+            const double zwoc = reg.zwoc ();
+            const double zgoc = reg.zgoc ();
+
+            if (! ((zgoc > z0) || (z0 > zwoc))) {
+                // Datum depth in oil zone  (zgoc <= z0 <= zwoc)
+                Details::equilibrateOWG(G, reg, grav, span, cells, press);
+            } else {
+                OPM_THROW(std::runtime_error, "Cannot initialise: the datum depth must be in the oil zone.");
+            }
+
+            return press;
+        }
+
+
+
+
+        template <class Region, class CellRange>
+        std::vector< std::vector<double> >
+        phaseSaturations(const Region&           reg,
+                         const CellRange&        cells,
+                         const BlackoilPropertiesInterface& props,
+                         std::vector< std::vector<double> >& phase_pressures)
+        {
+            const double z0   = reg.datum();
+            const double zwoc = reg.zwoc ();
+            const double zgoc = reg.zgoc ();
+            if ((zgoc > z0) || (z0 > zwoc)) {
+                OPM_THROW(std::runtime_error, "Cannot initialise: the datum depth must be in the oil zone.");
+            }
+            if (!reg.phaseUsage().phase_used[BlackoilPhases::Liquid]) {
+                OPM_THROW(std::runtime_error, "Cannot initialise: not handling water-gas cases.");
+            }
+
+            std::vector< std::vector<double> > phase_saturations = phase_pressures; // Just to get the right size.
+            double smin[BlackoilPhases::MaxNumPhases] = { 0.0 };
+            double smax[BlackoilPhases::MaxNumPhases] = { 0.0 };
+
+            const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
+            const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
+            const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+            const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua];
+            const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+            std::vector<double>::size_type local_index = 0;
+            for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) {
+                const int cell = *ci;
+                props.satRange(1, &cell, smin, smax);
+                // Find saturations from pressure differences by
+                // inverting capillary pressure functions.
+                double sw = 0.0;
+                if (water) {
+                    const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
+                    sw = satFromPc(props, waterpos, cell, pcov);
+                    phase_saturations[waterpos][local_index] = sw;
+                }
+                double sg = 0.0;
+                if (gas) {
+                    // Note that pcog is defined to be (pg - po), not (po - pg).
+                    const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index];
+                    const double increasing = true; // pcog(sg) expected to be increasing function
+                    sg = satFromPc(props, gaspos, cell, pcog, increasing);
+                    phase_saturations[gaspos][local_index] = sg;
+                }
+                bool overlap = false;
+                if (gas && water && (sg + sw > 1.0)) {
+                    // Overlapping gas-oil and oil-water transition
+                    // zones can lead to unphysical saturations when
+                    // treated as above. Must recalculate using gas-water
+                    // capillary pressure.
+                    const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index];
+                    sw = satFromSumOfPcs(props, waterpos, gaspos, cell, pcgw);
+                    sg = 1.0 - sw;
+                    phase_saturations[waterpos][local_index] = sw;
+                    phase_saturations[gaspos][local_index] = sg;
+                    overlap = true;
+                }
+                phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
+                
+                // Adjust phase pressures for max and min saturation ...
+                double pc[BlackoilPhases::MaxNumPhases];
+                double sat[BlackoilPhases::MaxNumPhases];
+                if (sw > smax[waterpos]-1.0e-6) {
+                    sat[waterpos] = smax[waterpos];
+                    props.capPress(1, sat, &cell, pc, 0);                   
+                    phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
+                } else if (overlap || sg > smax[gaspos]-1.0e-6) {
+                    sat[gaspos] = smax[gaspos];
+                    props.capPress(1, sat, &cell, pc, 0);                   
+                    phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
+                }
+                if (sg < smin[gaspos]+1.0e-6) {
+                    sat[gaspos] = smin[gaspos];
+                    props.capPress(1, sat, &cell, pc, 0);
+                    phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
+                }
+                if (sw < smin[waterpos]+1.0e-6) {
+                    sat[waterpos] = smin[waterpos];
+                    props.capPress(1, sat, &cell, pc, 0);
+                    phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
+                }
+            }
+            return phase_saturations;
+        }
+
+
+        /**
+         * Compute initial Rs values.
+         *
+         * \tparam CellRangeType Type of cell range that demarcates the
+         *                cells pertaining to the current
+         *                equilibration region.  Must implement
+         *                methods begin() and end() to bound the range
+         *                as well as provide an inner type,
+         *                const_iterator, to traverse the range.
+         *
+         * \param[in] grid            Grid.
+         * \param[in] cells           Range that spans the cells of the current
+         *                            equilibration region.
+         * \param[in] oil_pressure    Oil pressure for each cell in range.
+         * \param[in] rs_func         Rs as function of pressure and depth.
+         * \return                    Rs values, one for each cell in the 'cells' range.
+         */
+        template <class CellRangeType>
+        std::vector<double> computeRs(const UnstructuredGrid& grid,
+                                      const CellRangeType& cells,
+                                      const std::vector<double> oil_pressure,
+                                      const Miscibility::RsFunction& rs_func,
+                                      const std::vector<double> gas_saturation)
+        {
+            assert(grid.dimensions == 3);
+            std::vector<double> rs(cells.size());
+            int count = 0;
+            for (auto it = cells.begin(); it != cells.end(); ++it, ++count) {
+                const double depth = grid.cell_centroids[3*(*it) + 2];
+                rs[count] = rs_func(depth, oil_pressure[count], gas_saturation[count]);
+            }
+            return rs;
+        }
+
+    } // namespace Equil
+
+
+    namespace
+    {
+        /// Convert saturations from a vector of individual phase saturation vectors
+        /// to an interleaved format where all values for a given cell come before all
+        /// values for the next cell, all in a single vector.
+        std::vector<double> convertSats(const std::vector< std::vector<double> >& sat)
+        {
+            const int np = sat.size();
+            const int nc = sat[0].size();
+            std::vector<double> s(np * nc);
+            for (int c = 0; c < nc; ++c) {
+                for (int p = 0; p < np; ++p) {
+                    s[np*c + p] = sat[p][c];
+                }
+            }
+            return s;
+        }
+    }
+
+
+    /**
+     * Compute initial state by an equilibration procedure.
+     *
+     * The following state fields are modified:
+     *   pressure(),
+     *   saturation(),
+     *   surfacevol(),
+     *   gasoilratio(),
+     *   rv().
+     *
+     * \param[in] grid     Grid.
+     * \param[in] props    Property object, pvt and capillary properties are used.
+     * \param[in] deck     Simulation deck, used to obtain EQUIL and related data.
+     * \param[in] gravity  Acceleration of gravity, assumed to be in Z direction.
+     */
+    void initStateEquil(const UnstructuredGrid& grid,
+                        const BlackoilPropertiesInterface& props,
+                        const EclipseGridParser& deck,
+                        const double gravity,
+                        BlackoilState& state)
+    {
+        typedef Equil::DeckDependent::InitialStateComputer<EclipseGridParser> ISC;
+        ISC isc(props, deck, grid, gravity);
+        const auto pu = props.phaseUsage();
+        const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
+            ? pu.phase_pos[BlackoilPhases::Liquid]
+            : pu.phase_pos[BlackoilPhases::Aqua];
+        state.pressure() = isc.press()[ref_phase];
+        state.saturation() = convertSats(isc.saturation());
+        state.gasoilratio() = isc.rs();
+        state.rv() = isc.rv();
+        // TODO: state.surfacevol() must be computed from s, rs, rv.
+    }
+    
+    
+    void initStateEquil(const UnstructuredGrid& grid,
+                        const BlackoilPropertiesInterface& props,
+                        const Opm::DeckConstPtr newParserDeck,
+                        const double gravity,
+                        BlackoilState& state)
+    {
+        typedef Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> ISC;
+        ISC isc(props, newParserDeck, grid, gravity);
+        const auto pu = props.phaseUsage();
+        const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
+            ? pu.phase_pos[BlackoilPhases::Liquid]
+            : pu.phase_pos[BlackoilPhases::Aqua];
+        state.pressure() = isc.press()[ref_phase];
+        state.saturation() = convertSats(isc.saturation());
+        state.gasoilratio() = isc.rs();
+        state.rv() = isc.rv();
+        // TODO: state.surfacevol() must be computed from s, rs, rv.
+    }
+
+
+
+} // namespace Opm
+
+#endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED

tests/capillary.DATA

@@ -0,0 +1,38 @@
+WATER
+OIL
+GAS
+
+TABDIMS
+  1    1   40   20    1   20  /
+  
+EQLDIMS
+-- NTEQUL
+     1 /
+
+PVDO
+100 1.0 1.0
+200 0.9 1.0
+/
+
+PVDG
+100 0.010 0.1
+200 0.005 0.2
+/
+
+SWOF
+0.2 0 1 0.4
+1   1 0 0.1
+/
+
+SGOF
+0   0 1 0.2
+0.8 1 0 0.5
+/
+
+DENSITY
+700 1000 1
+/
+
+EQUIL
+50 150 50 0.25 20 0.35 1* 1* 0
+/

tests/capillary_overlap.DATA

@@ -0,0 +1,128 @@
+NOECHO
+
+RUNSPEC   ======
+
+WATER
+OIL
+GAS
+
+TABDIMS
+  1    1   40   20    1   20  /
+
+DIMENS
+1 1 20
+/
+
+WELLDIMS
+   30   10    2   30 /
+
+START
+   1 'JAN' 1990  /
+
+NSTACK
+   25 /
+
+EQLDIMS
+-- NTEQUL
+     1 / 
+     
+
+FMTOUT
+FMTIN
+
+GRID      ======
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+20*5.0
+/
+
+
+PORO
+20*0.2
+/
+
+
+PERMZ
+  20*1.0
+/
+
+PERMY
+20*100.0
+/
+
+PERMX
+20*100.0
+/
+
+BOX
+ 1 1 1 1 1 1 /
+
+TOPS
+0.0
+/
+
+PROPS     ======
+
+
+PVDO
+100 1.0 1.0
+200 0.9 1.0
+/
+
+PVDG
+100 0.010 0.1
+200 0.005 0.2
+/
+
+SWOF
+0.2 0 1 0.9
+1   1 0 0.1
+/
+
+SGOF
+0   0 1 0.2
+0.8 1 0 0.5
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+
+DENSITY
+700 1000 1
+/
+
+SOLUTION  ======
+
+EQUIL
+45 150 50 0.25 45 0.35 1* 1* 0
+/
+
+RPTSOL
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RESTART=2' /
+
+SUMMARY   ======
+RUNSUM
+
+SEPARATE
+
+SCHEDULE  ======
+
+RPTSCHED
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RESTART=3' 'NEWTON=2' /
+
+
+END

tests/deadfluids.DATA

@@ -0,0 +1,38 @@
+WATER
+OIL
+GAS
+
+TABDIMS
+  1    1   40   20    1   20  /
+  
+EQLDIMS
+-- NTEQUL
+     1 /
+
+PVDO
+100 1.0 1.0
+200 0.5 1.0
+/
+
+PVDG
+100 0.05 0.1
+200 0.02 0.2
+/
+
+SWOF
+0 0 1 0
+1 1 0 0
+/
+
+SGOF
+0 0 1 0
+1 1 0 0
+/
+
+DENSITY
+700 1000 10
+/
+
+EQUIL
+5 150 5 0 2 0 1* 1* 0
+/

tests/equil_livegas.DATA

@@ -0,0 +1,131 @@
+NOECHO
+
+RUNSPEC   ======
+
+WATER
+OIL
+GAS
+VAPOIL
+
+TABDIMS
+  1    1   40   20    1   20  /
+
+DIMENS
+1 1 20
+/
+
+WELLDIMS
+   30   10    2   30 /
+
+START
+   1 'JAN' 1990  /
+
+NSTACK
+   25 /
+
+EQLDIMS
+-- NTEQUL
+     1 / 
+     
+
+FMTOUT
+FMTIN
+
+GRID      ======
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+20*5.0
+/
+
+
+PORO
+20*0.2
+/
+
+
+PERMZ
+  20*1.0
+/
+
+PERMY
+20*100.0
+/
+
+PERMX
+20*100.0
+/
+
+BOX
+ 1 1 1 1 1 1 /
+
+TOPS
+0.0
+/
+
+PROPS     ======
+
+PVDO
+100 1.0 1.0
+200 0.9 1.0
+/
+
+PVTG
+--  Pg     Rv        Bg       Vg
+   100   0.0001       0.010      0.1
+         0.0          0.0104     0.1 /
+   200   0.0004       0.005      0.2
+         0.0          0.0054     0.2 /
+/
+
+SWOF
+0.2 0 1 0.9
+1   1 0 0.1
+/
+
+SGOF
+0   0 1 0.2
+0.8 1 0 0.5
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+
+DENSITY
+700 1000 1
+/
+
+SOLUTION  ======
+
+EQUIL
+45 150 50 0.25 45 0.35 1* 1* 0
+/
+
+RPTSOL
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
+
+SUMMARY   ======
+RUNSUM
+
+SEPARATE
+
+SCHEDULE  ======
+
+RPTSCHED
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
+
+
+END

tests/equil_liveoil.DATA

@@ -0,0 +1,140 @@
+NOECHO
+
+RUNSPEC   ======
+
+WATER
+OIL
+GAS
+DISGAS
+
+TABDIMS
+  1    1   40   20    1   20  /
+
+DIMENS
+1 1 20
+/
+
+WELLDIMS
+   30   10    2   30 /
+
+START
+   1 'JAN' 1990  /
+
+NSTACK
+   25 /
+
+EQLDIMS
+-- NTEQUL
+     1 / 
+     
+
+FMTOUT
+FMTIN
+
+GRID      ======
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+20*5.0
+/
+
+
+PORO
+20*0.2
+/
+
+
+PERMZ
+  20*1.0
+/
+
+PERMY
+20*100.0
+/
+
+PERMX
+20*100.0
+/
+
+BOX
+ 1 1 1 1 1 1 /
+
+TOPS
+0.0
+/
+
+PROPS     ======
+
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+         0          1.    1.0000     1.20  /
+        20         40.    1.0120     1.17  /
+        40         80.    1.0255     1.14  /
+        60        120.    1.0380     1.11  /
+        80        160.    1.0510     1.08  /
+       100        200.    1.0630     1.06  /
+       120        240.    1.0750     1.03  /
+       140        280.    1.0870     1.00  /
+       160        320.    1.0985      .98  /
+       180        360.    1.1100      .95  /
+       200        400.    1.1200      .94
+                  500.    1.1189      .94  /
+ /
+
+PVDG
+100 0.010 0.1
+200 0.005 0.2
+/
+
+SWOF
+0.2 0 1 0.9
+1   1 0 0.1
+/
+
+SGOF
+0   0 1 0.2
+0.8 1 0 0.5
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+
+DENSITY
+700 1000 1
+/
+
+SOLUTION  ======
+
+EQUIL
+45 150 50 0.25 45 0.35 1* 1* 0
+/
+
+RPTSOL
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
+
+SUMMARY   ======
+RUNSUM
+
+SEPARATE
+
+SCHEDULE  ======
+
+RPTSCHED
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
+
+
+END

tests/equil_liveoil_grid.DATA

@@ -0,0 +1,65 @@
+WATER
+OIL
+GAS
+DISGAS
+
+DIMENS
+1 1 20
+/
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+20*5.0
+/
+
+TOPS
+4*0.0
+/
+
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+         0          1.    1.0000     1.20  /
+       100         40.    1.0120     1.17  /
+       200         80.    1.0255     1.14  /
+       300        120.    1.0380     1.11  /
+       400        160.    1.0510     1.08  /
+       500        200.    1.0630     1.06  /
+       600        240.    1.0750     1.03  /
+       700        280.    1.0870     1.00  /
+       800        320.    1.0985      .98  /
+       900        360.    1.1100      .95  /
+      1000        400.    1.1200      .94
+                  500.    1.1189      .94  /
+ /
+
+
+PVDG
+100 0.010 0.1
+200 0.005 0.2
+/
+
+SWOF
+0.2 0 1 0.9
+1   1 0 0.1
+/
+
+SGOF
+0   0 1 0.2
+0.8 1 0 0.5
+/
+
+DENSITY
+700 1000 1
+/
+
+EQUIL
+45 150 50 0.25 45 0.35
+/

tests/equil_rsvd_and_rvvd.DATA

@@ -0,0 +1,151 @@
+NOECHO
+
+RUNSPEC   ======
+
+WATER
+OIL
+GAS
+DISGAS
+VAPOIL
+
+TABDIMS
+  1    1   40   20    1   20  /
+
+DIMENS
+1 1 20
+/
+
+WELLDIMS
+   30   10    2   30 /
+
+START
+   1 'JAN' 1990  /
+
+NSTACK
+   25 /
+
+EQLDIMS
+-- NTEQUL
+     1 / 
+     
+
+FMTOUT
+FMTIN
+
+GRID      ======
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+20*5.0
+/
+
+
+PORO
+20*0.2
+/
+
+
+PERMZ
+  20*1.0
+/
+
+PERMY
+20*100.0
+/
+
+PERMX
+20*100.0
+/
+
+BOX
+ 1 1 1 1 1 1 /
+
+TOPS
+0.0
+/
+
+PROPS     ======
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+         0          1.    1.0000     1.20  /
+        20         40.    1.0120     1.17  /
+        40         80.    1.0255     1.14  /
+        60        120.    1.0380     1.11  /
+        80        160.    1.0510     1.08  /
+       100        200.    1.0630     1.06  /
+       120        240.    1.0750     1.03  /
+       140        280.    1.0870     1.00  /
+       160        320.    1.0985      .98  /
+       180        360.    1.1100      .95  /
+       200        400.    1.1200      .94
+                  500.    1.1189      .94  /
+/
+
+PVTG
+--  Pg     Rv        Bg       Vg
+   100   0.0001       0.010      0.1
+         0.0          0.0104     0.1 /
+   200   0.0004       0.005      0.2
+         0.0          0.0054     0.2 /
+/
+
+SWOF
+0.2 0 1 0.9
+1   1 0 0.1
+/
+
+SGOF
+0   0 1 0.2
+0.8 1 0 0.5
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+
+DENSITY
+700 1000 1
+/
+
+SOLUTION  ======
+
+EQUIL
+45 150 50 0.25 45 0.35 1 1 0
+/
+
+RSVD
+ 0  0.0
+ 100 100. /
+ 
+RVVD
+   0.  0.
+ 100.  0.0001 /
+
+RPTSOL
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
+
+SUMMARY   ======
+RUNSUM
+
+SEPARATE
+
+SCHEDULE  ======
+
+RPTSCHED
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
+
+
+END

tests/test_equil.cpp

@@ -0,0 +1,764 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+*/
+
+#include "config.h"
+
+/* --- Boost.Test boilerplate --- */
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE  // Suppress own messages when throw()ing
+
+#define BOOST_TEST_MODULE UnitsTest
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+/* --- our own headers --- */
+
+#include <opm/core/simulator/initStateEquil.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/cart_grid.h>
+#include <opm/core/grid/GridManager.hpp>
+
+#include <opm/core/props/BlackoilPropertiesBasic.hpp>
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
+#include <opm/core/pressure/msmfem/partition.h>
+
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/Units.hpp>
+
+#include <array>
+#include <iostream>
+#include <limits>
+#include <memory>
+#include <numeric>
+#include <sstream>
+#include <string>
+#include <vector>
+
+BOOST_AUTO_TEST_SUITE ()
+
+BOOST_AUTO_TEST_CASE (PhasePressure)
+{
+    typedef std::vector<double> PVal;
+    typedef std::vector<PVal> PPress;
+
+    std::shared_ptr<UnstructuredGrid>
+        G(create_grid_cart3d(10, 1, 10), destroy_grid);
+
+    Opm::parameter::ParameterGroup param;
+    {
+        using Opm::unit::kilogram;
+        using Opm::unit::meter;
+        using Opm::unit::cubic;
+
+        std::stringstream dens; dens << 700*kilogram/cubic(meter);
+        param.insertParameter("rho2", dens.str());
+    }
+
+    typedef Opm::BlackoilPropertiesBasic Props;
+    Props props(param, G->dimensions, G->number_of_cells);
+
+    typedef Opm::Equil::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
+    RhoCalc calc(props, 0);
+
+    Opm::Equil::EquilRecord record =
+        {
+            { 0 , 1e5 } , // Datum depth, pressure
+            { 5 , 0   } , // Zwoc       , Pcow_woc
+            { 0 , 0   }   // Zgoc       , Pcgo_goc
+        };
+
+    Opm::Equil::EquilReg<RhoCalc>
+        region(record, calc,
+               std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+               std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+               props.phaseUsage());
+
+    std::vector<int> cells(G->number_of_cells);
+    std::iota(cells.begin(), cells.end(), 0);
+
+    const double grav   = 10;
+    const PPress ppress = Opm::Equil::phasePressures(*G, region, cells, grav);
+
+    const int first = 0, last = G->number_of_cells - 1;
+    const double reltol = 1.0e-8;
+    BOOST_CHECK_CLOSE(ppress[0][first] ,  90e3   , reltol);
+    BOOST_CHECK_CLOSE(ppress[0][last ] , 180e3   , reltol);
+    BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
+    BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
+}
+
+
+
+
+BOOST_AUTO_TEST_CASE (CellSubset)
+{
+    typedef std::vector<double> PVal;
+    typedef std::vector<PVal> PPress;
+
+    std::shared_ptr<UnstructuredGrid>
+        G(create_grid_cart3d(10, 1, 10), destroy_grid);
+
+    Opm::parameter::ParameterGroup param;
+    {
+        using Opm::unit::kilogram;
+        using Opm::unit::meter;
+        using Opm::unit::cubic;
+
+        std::stringstream dens; dens << 700*kilogram/cubic(meter);
+        param.insertParameter("rho2", dens.str());
+    }
+
+    typedef Opm::BlackoilPropertiesBasic Props;
+    Props props(param, G->dimensions, G->number_of_cells);
+
+    typedef Opm::Equil::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
+    RhoCalc calc(props, 0);
+
+    Opm::Equil::EquilRecord record[] =
+        {
+            {
+                { 0   ,  1e5     } , // Datum depth, pressure
+                { 2.5 , -0.075e5 } , // Zwoc       , Pcow_woc
+                { 0   ,  0       }   // Zgoc       , Pcgo_goc
+            }
+            ,
+            {
+                { 5   ,  1.35e5  } , // Datum depth, pressure
+                { 7.5 , -0.225e5 } , // Zwoc       , Pcow_woc
+                { 5   ,  0       }   // Zgoc       , Pcgo_goc
+            }
+        };
+
+    Opm::Equil::EquilReg<RhoCalc> region[] =
+        {
+            Opm::Equil::EquilReg<RhoCalc>(record[0], calc,
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          props.phaseUsage())
+            ,
+            Opm::Equil::EquilReg<RhoCalc>(record[0], calc,
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          props.phaseUsage())
+            ,
+            Opm::Equil::EquilReg<RhoCalc>(record[1], calc,
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          props.phaseUsage())
+            ,
+            Opm::Equil::EquilReg<RhoCalc>(record[1], calc,
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          props.phaseUsage())
+        };
+
+    const int cdim[] = { 2, 1, 2 };
+    int ncoarse = cdim[0];
+    for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; }
+
+    std::vector< std::vector<int> > cells(ncoarse);
+    for (int c = 0; c < G->number_of_cells; ++c) {
+        int ci = c;
+        const int i = ci % G->cartdims[0];  ci /= G->cartdims[0];
+        const int j = ci % G->cartdims[1];
+        const int k = ci / G->cartdims[1];
+
+        const int ic = (i / (G->cartdims[0] / cdim[0]));
+        const int jc = (j / (G->cartdims[1] / cdim[1]));
+        const int kc = (k / (G->cartdims[2] / cdim[2]));
+        const int ix = ic + cdim[0]*(jc + cdim[1]*kc);
+
+        assert ((0 <= ix) && (ix < ncoarse));
+        cells[ix].push_back(c);
+    }
+
+    PPress ppress(2, PVal(G->number_of_cells, 0));
+    for (std::vector< std::vector<int> >::const_iterator
+             r = cells.begin(), e = cells.end();
+         r != e; ++r)
+    {
+        const int    rno  = int(r - cells.begin());
+        const double grav = 10;
+        const PPress p    =
+            Opm::Equil::phasePressures(*G, region[rno], *r, grav);
+
+        PVal::size_type i = 0;
+        for (std::vector<int>::const_iterator
+                 c = r->begin(), ce = r->end();
+             c != ce; ++c, ++i)
+        {
+            assert (i < p[0].size());
+
+            ppress[0][*c] = p[0][i];
+            ppress[1][*c] = p[1][i];
+        }
+    }
+
+    const int first = 0, last = G->number_of_cells - 1;
+    const double reltol = 1.0e-8;
+    BOOST_CHECK_CLOSE(ppress[0][first] , 105e3   , reltol);
+    BOOST_CHECK_CLOSE(ppress[0][last ] , 195e3   , reltol);
+    BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
+    BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
+}
+
+
+
+
+BOOST_AUTO_TEST_CASE (RegMapping)
+{
+    typedef std::vector<double> PVal;
+    typedef std::vector<PVal> PPress;
+
+    std::shared_ptr<UnstructuredGrid>
+        G(create_grid_cart3d(10, 1, 10), destroy_grid);
+
+    Opm::parameter::ParameterGroup param;
+    {
+        using Opm::unit::kilogram;
+        using Opm::unit::meter;
+        using Opm::unit::cubic;
+
+        std::stringstream dens; dens << 700*kilogram/cubic(meter);
+        param.insertParameter("rho2", dens.str());
+    }
+
+    typedef Opm::BlackoilPropertiesBasic Props;
+    Props props(param, G->dimensions, G->number_of_cells);
+
+    typedef Opm::Equil::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
+    RhoCalc calc(props, 0);
+
+    Opm::Equil::EquilRecord record[] =
+        {
+            {
+                { 0   ,  1e5     } , // Datum depth, pressure
+                { 2.5 , -0.075e5 } , // Zwoc       , Pcow_woc
+                { 0   ,  0       }   // Zgoc       , Pcgo_goc
+            }
+            ,
+            {
+                { 5   ,  1.35e5  } , // Datum depth, pressure
+                { 7.5 , -0.225e5 } , // Zwoc       , Pcow_woc
+                { 5   ,  0       }   // Zgoc       , Pcgo_goc
+            }
+        };
+
+    Opm::Equil::EquilReg<RhoCalc> region[] =
+        {
+            Opm::Equil::EquilReg<RhoCalc>(record[0], calc,
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          props.phaseUsage())
+            ,
+            Opm::Equil::EquilReg<RhoCalc>(record[0], calc,
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          props.phaseUsage())
+            ,
+            Opm::Equil::EquilReg<RhoCalc>(record[1], calc,
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          props.phaseUsage())
+            ,
+            Opm::Equil::EquilReg<RhoCalc>(record[1], calc,
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
+                                          props.phaseUsage())
+        };
+
+    std::vector<int> eqlnum(G->number_of_cells);
+    {
+        std::vector<int> cells(G->number_of_cells);
+        std::iota(cells.begin(), cells.end(), 0);
+
+        const int cdim[] = { 2, 1, 2 };
+        int ncoarse = cdim[0];
+        for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; }
+
+        partition_unif_idx(G->dimensions, G->number_of_cells,
+                           G->cartdims, cdim,
+                           &cells[0], &eqlnum[0]);
+    }
+    Opm::RegionMapping<> eqlmap(eqlnum);
+
+    PPress ppress(2, PVal(G->number_of_cells, 0));
+    for (int r = 0, e = eqlmap.numRegions(); r != e; ++r)
+    {
+        const Opm::RegionMapping<>::CellRange&
+            rng = eqlmap.cells(r);
+
+        const int    rno  = r;
+        const double grav = 10;
+        const PPress p    =
+            Opm::Equil::phasePressures(*G, region[rno], rng, grav);
+
+        PVal::size_type i = 0;
+        for (Opm::RegionMapping<>::CellRange::const_iterator
+                 c = rng.begin(), ce = rng.end();
+             c != ce; ++c, ++i)
+        {
+            assert (i < p[0].size());
+
+            ppress[0][*c] = p[0][i];
+            ppress[1][*c] = p[1][i];
+        }
+    }
+
+    const int first = 0, last = G->number_of_cells - 1;
+    const double reltol = 1.0e-8;
+    BOOST_CHECK_CLOSE(ppress[0][first] , 105e3   , reltol);
+    BOOST_CHECK_CLOSE(ppress[0][last ] , 195e3   , reltol);
+    BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
+    BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
+}
+
+
+
+BOOST_AUTO_TEST_CASE (DeckAllDead)
+{
+    std::shared_ptr<UnstructuredGrid>
+        grid(create_grid_cart3d(1, 1, 10), destroy_grid);
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("deadfluids.DATA");
+    Opm::BlackoilPropertiesFromDeck props(deck, *grid, false);
+    Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, *grid, 10.0);
+    const auto& pressures = comp.press();
+    BOOST_REQUIRE(pressures.size() == 3);
+    BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells);
+
+    const int first = 0, last = grid->number_of_cells - 1;
+    // The relative tolerance is too loose to be very useful,
+    // but the answer we are checking is the result of an ODE
+    // solver, and it is unclear if we should check it against
+    // the true answer or something else.
+    const double reltol = 1.0e-3;
+    BOOST_CHECK_CLOSE(pressures[0][first] , 1.496329839e7   , reltol);
+    BOOST_CHECK_CLOSE(pressures[0][last ] , 1.50473245e7   , reltol);
+    BOOST_CHECK_CLOSE(pressures[1][last] , 1.50473245e7   , reltol);
+}
+
+
+
+BOOST_AUTO_TEST_CASE (CapillaryInversion)
+{
+    // Test setup.
+    Opm::GridManager gm(1, 1, 40, 1.0, 1.0, 2.5);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("capillary.DATA");
+    Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
+
+    // Test the capillary inversion for oil-water.
+    const int cell = 0;
+    const double reltol = 1.0e-7;
+    {
+        const int phase = 0;
+        const bool increasing = false;
+        const std::vector<double> pc = { 10.0e5, 0.5e5, 0.4e5, 0.3e5, 0.2e5, 0.1e5, 0.099e5, 0.0e5, -10.0e5 };
+        const std::vector<double> s = { 0.2, 0.2, 0.2, 0.466666666666, 0.733333333333, 1.0, 1.0, 1.0, 1.0 };
+        BOOST_REQUIRE(pc.size() == s.size());
+        for (size_t i = 0; i < pc.size(); ++i) {
+            const double s_computed = Opm::Equil::satFromPc(props, phase, cell, pc[i], increasing);
+            BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
+        }
+    }
+
+    // Test the capillary inversion for gas-oil.
+    {
+        const int phase = 2;
+        const bool increasing = true;
+        const std::vector<double> pc = { 10.0e5, 0.6e5, 0.5e5, 0.4e5, 0.3e5, 0.2e5, 0.1e5, 0.0e5, -10.0e5 };
+        const std::vector<double> s = { 0.8, 0.8, 0.8, 0.533333333333, 0.266666666666, 0.0, 0.0, 0.0, 0.0 };
+        BOOST_REQUIRE(pc.size() == s.size());
+        for (size_t i = 0; i < pc.size(); ++i) {
+            const double s_computed = Opm::Equil::satFromPc(props, phase, cell, pc[i], increasing);
+            BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
+        }
+    }
+
+    // Test the capillary inversion for gas-water.
+    {
+        const int water = 0;
+        const int gas = 2;
+        const std::vector<double> pc = { 0.9e5, 0.8e5, 0.6e5, 0.4e5, 0.3e5 };
+        const std::vector<double> s = { 0.2, 0.333333333333, 0.6, 0.866666666666, 1.0 };
+        BOOST_REQUIRE(pc.size() == s.size());
+        for (size_t i = 0; i < pc.size(); ++i) {
+            const double s_computed = Opm::Equil::satFromSumOfPcs(props, water, gas, cell, pc[i]);
+            BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
+        }
+    }
+}
+
+
+
+BOOST_AUTO_TEST_CASE (DeckWithCapillary)
+{
+    Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("capillary.DATA");
+    Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
+
+    Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 10.0);
+    const auto& pressures = comp.press();
+    BOOST_REQUIRE(pressures.size() == 3);
+    BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
+
+    const int first = 0, last = grid.number_of_cells - 1;
+    // The relative tolerance is too loose to be very useful,
+    // but the answer we are checking is the result of an ODE
+    // solver, and it is unclear if we should check it against
+    // the true answer or something else.
+    const double reltol = 1.0e-6;
+    BOOST_CHECK_CLOSE(pressures[0][first] , 1.469769063e7   , reltol);
+    BOOST_CHECK_CLOSE(pressures[0][last ] , 1.545e7   , reltol);
+    BOOST_CHECK_CLOSE(pressures[1][last] , 1.546e7   , reltol);
+
+    const auto& sats = comp.saturation();
+    const std::vector<double> s[3]{
+        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.425893333333, 0.774026666666, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+        { 0, 0, 0, 0.00736, 0.792746666666, 0.8, 0.8, 0.8, 0.8, 0.574106666666, 0.225973333333, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+        { 0.8, 0.8, 0.8, 0.79264, 0.007253333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
+    };
+    for (int phase = 0; phase < 3; ++phase) {
+        BOOST_REQUIRE(sats[phase].size() == s[phase].size());
+        for (size_t i = 0; i < s[phase].size(); ++i) {
+            BOOST_CHECK_CLOSE(sats[phase][i], s[phase][i], reltol);
+        }
+    }
+}
+
+
+
+BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
+{
+    Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("capillary_overlap.DATA");
+    Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
+
+    Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
+    const auto& pressures = comp.press();
+    BOOST_REQUIRE(pressures.size() == 3);
+    BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
+
+    const int first = 0, last = grid.number_of_cells - 1;
+    // The relative tolerance is too loose to be very useful,
+    // but the answer we are checking is the result of an ODE
+    // solver, and it is unclear if we should check it against
+    // the true answer or something else.
+    const double reltol = 1.0e-6;
+    const double reltol_ecl = 1.0;
+    BOOST_CHECK_CLOSE(pressures[0][first], 1.48324e+07, reltol_ecl);  // eclipse 
+    BOOST_CHECK_CLOSE(pressures[0][last],  1.54801e+07, reltol_ecl);
+    BOOST_CHECK_CLOSE(pressures[1][first], 1.49224e+07, reltol_ecl);
+    BOOST_CHECK_CLOSE(pressures[1][last],  1.54901e+07, reltol_ecl);
+    
+    BOOST_CHECK_CLOSE(pressures[0][first] , 14832467.14, reltol); // opm
+    BOOST_CHECK_CLOSE(pressures[0][last ] , 15479883.47, reltol);
+    BOOST_CHECK_CLOSE(pressures[1][last ] , 15489883.47, reltol);
+
+    const auto& sats = comp.saturation();
+    // std::cout << "Saturations:\n";
+    // for (const auto& sat : sats) {
+    //     for (const double s : sat) {
+    //         std::cout << s << ' ';
+    //     }
+    //     std::cout << std::endl;
+    // }
+   
+    const std::vector<double> s_ecl[3]{// eclipse
+        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22874042, 0.53397995, 0.78454906,  0.91542006, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+        { 0,   0,   0,   0,   0,   0,   0,   0,          0,          0.20039,     0.08458,    0, 0, 0, 0, 0, 0, 0, 0, 0 },
+        { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77125955, 0.46602005, 0.015063271, 0,          0, 0, 0, 0, 0, 0, 0, 0, 0 }
+    };
+
+    const std::vector<double> s_opm[3]{ // opm 
+        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2289309090909091,  0.53406545454545451, 0.78458,              0.9154, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+        { 0,   0,   0,   0,   0,   0,   0,   0,                   0,                   0.2002466666666666,   0.0846, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+        { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77106909090909093, 0.46593454545454549, 0.015173333333333336, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0 }
+    };
+    for (int phase = 0; phase < 3; ++phase) {
+        BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size());
+        for (size_t i = 0; i < s_opm[phase].size(); ++i) {
+            //std::cout << std::setprecision(10) << sats[phase][i] << '\n';
+            BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
+            BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
+        }
+    }
+}
+
+
+
+BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
+{
+    Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("equil_liveoil.DATA");
+    Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
+
+    Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
+    const auto& pressures = comp.press();
+    BOOST_REQUIRE(pressures.size() == 3);
+    BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
+
+    const int first = 0, last = grid.number_of_cells - 1;
+    // The relative tolerance is too loose to be very useful,
+    // but the answer we are checking is the result of an ODE
+    // solver, and it is unclear if we should check it against
+    // the true answer or something else.
+    const double reltol = 1.0e-6;
+    const double reltol_ecl = 1.0;
+    BOOST_CHECK_CLOSE(pressures[0][first], 1.48324e+07, reltol_ecl);  // eclipse
+    BOOST_CHECK_CLOSE(pressures[0][last],  1.54801e+07, reltol_ecl);
+    BOOST_CHECK_CLOSE(pressures[1][first], 1.49224e+07, reltol_ecl);
+    BOOST_CHECK_CLOSE(pressures[1][last],  1.54901e+07, reltol_ecl);
+    
+    BOOST_CHECK_CLOSE(pressures[0][first], 1.483246714e7, reltol);  // opm
+    BOOST_CHECK_CLOSE(pressures[0][last],  1.547991652e7, reltol);
+    BOOST_CHECK_CLOSE(pressures[1][first], 1.492246714e7, reltol);
+    BOOST_CHECK_CLOSE(pressures[1][last],  1.548991652e7, reltol);
+
+    const auto& sats = comp.saturation();
+    // std::cout << "Saturations:\n";
+    // for (const auto& sat : sats) {
+    //     for (const double s : sat) {
+    //         std::cout << s << ' ';
+    //     }
+    //     std::cout << std::endl;
+    // }
+    const std::vector<double> s_ecl[3]{ // eclipse
+        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22898, 0.53422, 0.78470, 0.91531, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+        { 0,   0,   0,   0,   0,   0,   0,   0,       0,       0.20073, 0.08469, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+        { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77102, 0.46578, 0.01458, 0,       0, 0, 0, 0, 0, 0, 0, 0, 0 }
+    };
+    const std::vector<double> s_opm[3]{ // opm 
+        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2291709091, 0.5343054545, 0.78472,        0.91529, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+        { 0,   0,   0,   0,   0,   0,   0,   0,            0,            0.2005866667,   0.08471, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+        { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.7708290909, 0.4656945455, 0.01469333333,  0,       0, 0, 0, 0, 0, 0, 0, 0, 0 }
+    };
+    for (int phase = 0; phase < 3; ++phase) {
+        BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size());
+        for (size_t i = 0; i < s_opm[phase].size(); ++i) {
+            //std::cout << std::setprecision(10) << sats[phase][i] << '\n';
+            BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
+            BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
+        }
+        std::cout << std::endl;
+    }
+    
+    const auto& rs = comp.rs();
+    const std::vector<double> rs_opm {74.612335679539058, 74.649052116644228, 74.685786561426298, 74.722539022717172, // opm
+                                      74.759309509353145, 74.796098030174733, 74.8329045940269,   74.869729209758916, 
+                                      74.906571886224327, 75.090675116639048, 75.0,               75.0, 
+                                      75.0,               75.0,               75.0,               75.0, 
+                                      75.0,               75.0,               75.0,               75.0};
+    const std::vector<double> rs_ecl {74.612228, 74.648956, 74.685707, 74.722473,  // eclipse
+                                      74.759254, 74.796051, 74.832870, 74.875145,
+                                      74.969231, 75.090706, 75.000000, 75.000000,
+                                      75.000000, 75.000000, 75.000000, 75.000000,
+                                      75.000000, 75.000000, 75.000000, 75.000000}; 
+    for (size_t i = 0; i < rs_opm.size(); ++i) {
+        //std::cout << std::setprecision(10) << sats[phase][i] << '\n';
+        BOOST_CHECK_CLOSE(rs[i], rs_opm[i], reltol);
+        BOOST_CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl);
+    }
+}
+
+
+
+BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
+{
+    Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("equil_livegas.DATA");
+    Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
+
+    Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
+    const auto& pressures = comp.press();
+    BOOST_REQUIRE(pressures.size() == 3);
+    BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
+
+    const int first = 0, last = grid.number_of_cells - 1;
+    // The relative tolerance is too loose to be very useful,
+    // but the answer we are checking is the result of an ODE
+    // solver, and it is unclear if we should check it against
+    // the true answer or something else.
+    const double reltol = 1.0e-6;
+    const double reltol_ecl = 1.0;
+    BOOST_CHECK_CLOSE(pressures[0][first], 1.48215e+07, reltol_ecl);  // eclipse
+    BOOST_CHECK_CLOSE(pressures[0][last],  1.54801e+07, reltol_ecl);
+    BOOST_CHECK_CLOSE(pressures[1][first], 1.49115e+07, reltol_ecl);
+    BOOST_CHECK_CLOSE(pressures[1][last],  1.54901e+07, reltol_ecl);
+    
+    BOOST_CHECK_CLOSE(pressures[0][first], 1.482150311e7, reltol);  // opm
+    BOOST_CHECK_CLOSE(pressures[0][last],  1.547988347e7, reltol);
+    BOOST_CHECK_CLOSE(pressures[1][first], 1.491150311e7, reltol);
+    BOOST_CHECK_CLOSE(pressures[1][last],  1.548988347e7, reltol);
+
+    const auto& sats = comp.saturation();
+    // std::cout << "Saturations:\n";
+    // for (const auto& sat : sats) {
+    //     for (const double s : sat) {
+    //         std::cout << s << ' ';
+    //     }
+    //     std::cout << std::endl;
+    // }
+    const std::vector<double> s_ecl[3]{ // eclipse
+        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.24285614, 0.53869015, 0.78454906,  0.91542006, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+        { 0,   0,   0,   0,   0,   0,   0,   0,          0,          0.18311,     0.08458,    0, 0, 0, 0, 0, 0, 0, 0, 0 },
+        { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.75714386, 0.46130988, 0.032345835, 0,          0, 0, 0, 0, 0, 0, 0, 0, 0 }
+    };
+
+    const std::vector<double> s_opm[3]{ // opm 
+        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.24310545, 0.5388, 0.78458,    0.91540, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+        { 0,   0,   0,   0,   0,   0,   0,   0,          0,      0.18288667, 0.0846,  0, 0, 0, 0, 0, 0, 0, 0, 0 },
+        { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.75689455, 0.4612, 0.03253333, 0,       0, 0, 0, 0, 0, 0, 0, 0, 0 }
+    };
+    for (int phase = 0; phase < 3; ++phase) {
+        BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size());
+        for (size_t i = 0; i < s_opm[phase].size(); ++i) {
+            //std::cout << std::setprecision(10) << sats[phase][i] << '\n';
+            BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], 100.*reltol);
+            BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
+        }
+        std::cout << std::endl;
+    }
+    
+    const auto& rv = comp.rv();
+    const std::vector<double> rv_opm { // opm
+        2.4884509e-4, 2.4910378e-4, 2.4936267e-4, 2.4962174e-4,
+        2.4988100e-4, 2.5014044e-4, 2.5040008e-4, 2.5065990e-4, 
+        2.5091992e-4, 2.5118012e-4, 2.5223082e-4, 2.5105e-4, 
+        2.5105e-4,    2.5105e-4,    2.5105e-4,    2.5105e-4, 
+        2.5105e-4,    2.5105e-4,    2.5105e-4,    2.5105e-4};
+
+    const std::vector<double> rv_ecl {  // eclipse
+        0.24884584E-03,   0.24910446E-03,   0.24936325E-03,   0.24962222E-03,
+        0.24988138E-03,   0.25014076E-03,   0.25040031E-03,   0.25066003E-03,
+        0.25091995E-03,   0.25118008E-03,   0.25223137E-03,   0.25104999E-03,
+        0.25104999E-03,   0.25104999E-03,   0.25104999E-03,   0.25104999E-03,
+        0.25104999E-03,   0.25104999E-03,   0.25104999E-03,   0.25104999E-03};
+         
+    for (size_t i = 0; i < rv_opm.size(); ++i) {
+        //std::cout << std::setprecision(10) << sats[phase][i] << '\n';
+        BOOST_CHECK_CLOSE(rv[i], rv_opm[i], 100.*reltol);
+        BOOST_CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl);
+    }
+}
+
+BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
+{
+    Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("equil_rsvd_and_rvvd.DATA");
+    Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
+
+    Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
+    const auto& pressures = comp.press();
+    BOOST_REQUIRE(pressures.size() == 3);
+    BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
+
+    const int first = 0, last = grid.number_of_cells - 1;
+    // The relative tolerance is too loose to be very useful,
+    // but the answer we are checking is the result of an ODE
+    // solver, and it is unclear if we should check it against
+    // the true answer or something else.
+    const double reltol = 1.0e-6;
+    const double reltol_ecl = 1.0;
+    BOOST_CHECK_CLOSE(pressures[0][first], 1.48350e+07, reltol_ecl);  // eclipse
+    BOOST_CHECK_CLOSE(pressures[0][last],  1.54794e+07, reltol_ecl);
+    BOOST_CHECK_CLOSE(pressures[1][first], 1.49250e+07, reltol_ecl);
+    BOOST_CHECK_CLOSE(pressures[1][last],  1.54894e+07, reltol_ecl);
+    
+    BOOST_CHECK_CLOSE(pressures[0][first], 1.483499660e7, reltol);  // opm
+    BOOST_CHECK_CLOSE(pressures[0][last],  1.547924516e7, reltol);
+    BOOST_CHECK_CLOSE(pressures[1][first], 1.492499660e7, reltol);
+    BOOST_CHECK_CLOSE(pressures[1][last],  1.548924516e7, reltol);
+
+    const auto& sats = comp.saturation();
+    // std::cout << "Saturations:\n";
+    // for (const auto& sat : sats) {
+    //     for (const double s : sat) {
+    //         std::cout << s << ' ';
+    //     }
+    //     std::cout << std::endl;
+    // }
+    const std::vector<double> s_ecl[3]{ // eclipse
+        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22206347, 0.52871972, 0.78150368,  0.91819441,  1, 1, 1, 1, 1, 1, 1, 1, 1 },
+        { 0,   0,   0,   0,   0,   0,   0,   0,          0,          0.19656529,  0.081805572, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+        { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77793652, 0.47128031, 0.021931054, 0,           0, 0, 0, 0, 0, 0, 0, 0, 0 }
+    };
+
+    const std::vector<double> s_opm[3]{ // opm 
+        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22232000, 0.52882909, 0.78153000,    0.91817000, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+        { 0,   0,   0,   0,   0,   0,   0,   0,          0,          0.19636333,    0.08183000, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+        { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77768000, 0.47117091, 0.02210667,    0,          0, 0, 0, 0, 0, 0, 0, 0, 0 }
+    };
+    
+    for (int phase = 0; phase < 3; ++phase) {
+        BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size());
+        for (size_t i = 0; i < s_opm[phase].size(); ++i) {
+            //std::cout << std::setprecision(10) << sats[phase][i] << '\n';
+            BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], 100.*reltol);
+            BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
+        }
+        std::cout << std::endl;
+    }
+    
+    const auto& rs = comp.rs();
+    const std::vector<double> rs_opm { // opm
+        74.624983020822540, 74.659590408801634, 74.694380353364522, 74.729353362649505,
+        74.764509945812975, 74.799850613032362, 74.835375875509555, 74.87108624547416, 
+        74.906982236186707, 75.088917653469309, 52.5,               57.5, 
+        62.5,               67.5,               72.5,               76.45954840804761, 
+        76.70621044909619,  76.952877357524045, 77.199549133522638, 77.446225777283587};
+
+    const std::vector<double> rs_ecl {  // eclipse
+        74.625114, 74.659706, 74.694481, 74.729439,
+        74.764580, 74.799904, 74.835419, 74.875252,
+        74.968628, 75.088951, 52.500000, 57.500000,
+        62.500000, 67.500000, 72.500000, 76.168388,
+        76.349953, 76.531532, 76.713142, 76.894775,};
+    
+    const auto& rv = comp.rv();
+    const std::vector<double> rv_opm { // opm
+        2.50e-6, 7.50e-6,       1.25e-5,       1.75e-5,
+        2.25e-5, 2.75e-5,       3.25e-5,       3.75e-5, 
+        4.25e-5, 2.51158386e-4, 2.52203372e-4, 5.75e-5, 
+        6.25e-5, 6.75e-5,       7.25e-5,       7.75e-5, 
+        8.25e-5, 8.75e-5,       9.25e-5,       9.75e-5};
+
+    const std::vector<double> rv_ecl {  // eclipse
+        0.24999999E-05, 0.74999998E-05, 0.12500000E-04, 0.17500000E-04,
+        0.22500000E-04, 0.27500000E-04, 0.32500000E-04, 0.37500002E-04,
+        0.42500000E-04, 0.25115837E-03, 0.25220393E-03, 0.57500001E-04,
+        0.62500003E-04, 0.67499997E-04, 0.72499999E-04, 0.77500001E-04,
+        0.82500002E-04, 0.87499997E-04, 0.92499999E-04, 0.97500000E-04};
+         
+    for (size_t i = 0; i < rv_opm.size(); ++i) {
+        //std::cout << std::setprecision(10) << sats[phase][i] << '\n';
+        BOOST_CHECK_CLOSE(rs[i], rs_opm[i], 100*reltol);
+        BOOST_CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl);
+        BOOST_CHECK_CLOSE(rv[i], rv_opm[i], 100.*reltol);
+        BOOST_CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl);
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()