diff --git a/examples/compute_initial_state.cpp b/examples/compute_initial_state.cpp
new file mode 100644
index 000000000..3193c9fac
--- /dev/null
+++ b/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 .
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+
+#include
+
+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& 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(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("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("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;
+}
diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp
new file mode 100644
index 000000000..d89dd8c3c
--- /dev/null
+++ b/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 .
+*/
+
+#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
+#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+
+/*
+---- synopsis of EquilibrationHelpers.hpp ----
+
+namespace Opm
+{
+ namespace Equil {
+
+ template
+ class DensityCalculator;
+
+ template <>
+ class DensityCalculator< BlackoilPropertiesInterface >;
+
+ namespace Miscibility {
+ class RsFunction;
+ class NoMixing;
+ class RsVD;
+ class RsSatAtContact;
+ }
+
+ struct EquilRecord;
+
+ template
+ 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 DensityCalculator;
+
+ /**
+ * Facility for calculating phase densities based on the
+ * BlackoilPropertiesInterface.
+ *
+ * Implements the crucial operator()(p,svol)
+ * 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
+ operator()(const double p,
+ const std::vector& z) const
+ {
+ const int np = props_.numPhases();
+ std::vector A(np * np, 0);
+
+ assert (z.size() == std::vector::size_type(np));
+
+ double* dAdp = 0;
+ props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
+
+ std::vector rho(np, 0.0);
+ props_.density(1, &A[0], &rho[0]);
+
+ return rho;
+ }
+
+ private:
+ const BlackoilPropertiesInterface& props_;
+ const std::vector 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& depth,
+ const std::vector& 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 depth_; /**< Depth nodes */
+ std::vector 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& depth,
+ const std::vector& 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 depth_; /**< Depth nodes */
+ std::vector 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
+ *
+ * std::vector
+ * operator()(const double press,
+ * const std::vector& svol )
+ *
+ * that calculates the phase densities of all phases in @c
+ * svol at fluid pressure @c press.
+ */
+ template
+ 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 rs,
+ std::shared_ptr 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 rs_; /**< RS calculator */
+ std::shared_ptr 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 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 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
diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp
new file mode 100644
index 000000000..c1f72b714
--- /dev/null
+++ b/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 .
+*/
+
+#ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
+#define OPM_INITSTATEEQUIL_HEADER_INCLUDED
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+/**
+ * \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
+ std::vector< std::vector >
+ 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
+ std::vector< std::vector >
+ phaseSaturations(const Region& reg,
+ const CellRange& cells,
+ const BlackoilPropertiesInterface& props,
+ const std::vector< std::vector >& 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
+ std::vector computeRs(const UnstructuredGrid& grid,
+ const CellRangeType& cells,
+ const std::vector oil_pressure,
+ const Miscibility::RsFunction& rs_func,
+ const std::vector gas_saturation);
+
+ namespace DeckDependent {
+
+ inline
+ std::vector
+ getEquil(const EclipseGridParser& deck)
+ {
+ if (deck.hasField("EQUIL")) {
+ const EQUIL& eql = deck.getEQUIL();
+
+ typedef std::vector::size_type sz_t;
+ const sz_t nrec = eql.equil.size();
+
+ std::vector 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
+ getEquil(const Opm::DeckConstPtr newParserDeck)
+ {
+ if (newParserDeck->hasKeyword("EQUIL")) {
+
+ Opm::EquilWrapper eql(newParserDeck->getKeyword("EQUIL"));
+
+ const int nrec = eql.numRegions();
+
+ std::vector 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
+ equilnum(const EclipseGridParser& deck,
+ const UnstructuredGrid& G )
+ {
+ std::vector eqlnum;
+ if (deck.hasField("EQLNUM")) {
+ const std::vector& e = deck.getIntegerValue("EQLNUM");
+ eqlnum.reserve(e.size());
+ std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
+ std::bind2nd(std::minus(), 1));
+ }
+ else {
+ // No explicit equilibration region.
+ // All cells in region zero.
+ eqlnum.assign(G.number_of_cells, 0);
+ }
+
+ return eqlnum;
+ }
+
+
+ inline
+ std::vector
+ equilnum(const Opm::DeckConstPtr newParserDeck,
+ const UnstructuredGrid& G )
+ {
+ std::vector eqlnum;
+ if (newParserDeck->hasKeyword("EQLNUM")) {
+ const std::vector& e = newParserDeck->getKeyword("EQLNUM")->getIntData();
+ eqlnum.reserve(e.size());
+ std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
+ std::bind2nd(std::minus(), 1));
+ }
+ else {
+ // No explicit equilibration region.
+ // All cells in region zero.
+ eqlnum.assign(G.number_of_cells, 0);
+ }
+
+ return eqlnum;
+ }
+
+
+ template
+ class InitialStateComputer;
+
+ template <>
+ class InitialStateComputer {
+ public:
+ InitialStateComputer(const BlackoilPropertiesInterface& props,
+ const EclipseGridParser& deck ,
+ const UnstructuredGrid& G ,
+ const double grav = unit::gravity)
+ : pp_(props.numPhases(),
+ std::vector(G.number_of_cells)),
+ sat_(props.numPhases(),
+ std::vector(G.number_of_cells)),
+ rs_(G.number_of_cells),
+ rv_(G.number_of_cells)
+ {
+ // Get the equilibration records.
+ const std::vector 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 depth; depth.push_back(0.0); depth.push_back(100.0);
+ std::vector rs; rs.push_back(0.0); rs.push_back(100.0);
+ rs_func_.push_back(std::make_shared(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(props, cell, p_contact));
+ }
+ }
+ } else {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ rs_func_.push_back(std::make_shared());
+ }
+ }
+
+ 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 depth; depth.push_back(0.0); depth.push_back(100.0);
+ std::vector rv; rv.push_back(0.0); rv.push_back(0.0001);
+ rv_func_.push_back(std::make_shared(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(props, cell, p_contact));
+ }
+ }
+ } else {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ rv_func_.push_back(std::make_shared());
+ }
+ }
+
+ // 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 Vec;
+ typedef std::vector 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 RhoCalc;
+ typedef EquilReg EqReg;
+
+ std::vector< std::shared_ptr > rs_func_;
+ std::vector< std::shared_ptr > rv_func_;
+
+ PVec pp_;
+ PVec sat_;
+ Vec rs_;
+ Vec rv_;
+
+ template
+ 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
+ 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 {
+ public:
+ InitialStateComputer(const BlackoilPropertiesInterface& props,
+ const Opm::DeckConstPtr newParserDeck,
+ const UnstructuredGrid& G ,
+ const double grav = unit::gravity)
+ : pp_(props.numPhases(),
+ std::vector(G.number_of_cells)),
+ sat_(props.numPhases(),
+ std::vector(G.number_of_cells)),
+ rs_(G.number_of_cells),
+ rv_(G.number_of_cells)
+ {
+ // Get the equilibration records.
+ const std::vector 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{"vd", "rs"},rec[i].live_oil_table_index-1);
+ std::vector vd(rsvd.getColumn("vd"));
+ std::vector rs(rsvd.getColumn("rs"));
+ rs_func_.push_back(std::make_shared(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(props, cell, p_contact));
+ }
+ }
+ } else {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ rs_func_.push_back(std::make_shared());
+ }
+ }
+
+ 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{"vd", "rv"},rec[i].wet_gas_table_index-1);
+ std::vector vd(rvvd.getColumn("vd"));
+ std::vector rv(rvvd.getColumn("rv"));
+ rv_func_.push_back(std::make_shared(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(props, cell, p_contact));
+ }
+ }
+ } else {
+ for (size_t i = 0; i < rec.size(); ++i) {
+ rv_func_.push_back(std::make_shared());
+ }
+ }
+
+ // 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 Vec;
+ typedef std::vector 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 RhoCalc;
+ typedef EquilReg EqReg;
+
+ std::vector< std::shared_ptr > rs_func_;
+ std::vector< std::shared_ptr > rv_func_;
+
+ PVec pp_;
+ PVec sat_;
+ Vec rs_;
+ Vec rv_;
+
+ template
+ 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
+ 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
+
+#endif // OPM_INITSTATEEQUIL_HEADER_INCLUDED
diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp
new file mode 100644
index 000000000..d227d54d1
--- /dev/null
+++ b/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 .
+*/
+
+#ifndef OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
+#define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+namespace Opm
+{
+ namespace Details {
+ template
+ class RK4IVP {
+ public:
+ RK4IVP(const RHS& f ,
+ const std::array& 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::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 span_;
+ std::vector y_;
+ std::vector f_;
+
+ double
+ stepsize() const { return (span_[1] - span_[0]) / N_; }
+ };
+
+ namespace PhasePressODE {
+ template
+ 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 svol_;
+ const int ix_;
+ const double g_;
+
+ double
+ density(const double press) const
+ {
+ const std::vector& rho = rho_(press, svol_);
+
+ return rho[ix_];
+ }
+ };
+
+ template
+ 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 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& rho = rho_(press, svol_);
+ return rho[oix_];
+ }
+ };
+
+ template
+ 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 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& 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
+ void
+ assign(const UnstructuredGrid& G ,
+ const std::array& f ,
+ const double split,
+ const CellRange& cells,
+ std::vector& p )
+ {
+ const int nd = G.dimensions;
+
+ enum { up = 0, down = 1 };
+
+ std::vector::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
+ void
+ water(const UnstructuredGrid& G ,
+ const Region& reg ,
+ const std::array& span ,
+ const double grav ,
+ const double po_woc,
+ const CellRange& cells ,
+ std::vector& press )
+ {
+ using PhasePressODE::Water;
+ typedef Water 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 up = {{ z0, span[0] }};
+ std::array down = {{ z0, span[1] }};
+
+ typedef Details::RK4IVP WPress;
+ std::array wpress = {
+ {
+ WPress(drho, up , p0, 100)
+ ,
+ WPress(drho, down, p0, 100)
+ }
+ };
+
+ assign(G, wpress, z0, cells, press);
+ }
+
+ template
+ void
+ oil(const UnstructuredGrid& G ,
+ const Region& reg ,
+ const std::array& span ,
+ const double grav ,
+ const CellRange& cells ,
+ std::vector& press ,
+ double& po_woc,
+ double& po_goc)
+ {
+ using PhasePressODE::Oil;
+ typedef Oil 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 up = {{ z0, span[0] }};
+ std::array down = {{ z0, span[1] }};
+
+ typedef Details::RK4IVP OPress;
+ std::array opress = {
+ {
+ OPress(drho, up , p0, 100)
+ ,
+ OPress(drho, down, p0, 100)
+ }
+ };
+
+ assign(G, opress, z0, cells, press);
+
+ po_woc = -std::numeric_limits::max();
+ po_goc = -std::numeric_limits::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
+ void
+ gas(const UnstructuredGrid& G ,
+ const Region& reg ,
+ const std::array& span ,
+ const double grav ,
+ const double po_goc,
+ const CellRange& cells ,
+ std::vector& press )
+ {
+ using PhasePressODE::Gas;
+ typedef Gas 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 up = {{ z0, span[0] }};
+ std::array down = {{ z0, span[1] }};
+
+ typedef Details::RK4IVP GPress;
+ std::array gpress = {
+ {
+ GPress(drho, up , p0, 100)
+ ,
+ GPress(drho, down, p0, 100)
+ }
+ };
+
+ assign(G, gpress, z0, cells, press);
+ }
+ } // namespace PhasePressure
+
+ template
+ void
+ equilibrateOWG(const UnstructuredGrid& G,
+ const Region& reg,
+ const double grav,
+ const std::array& span,
+ const CellRange& cells,
+ std::vector< std::vector >& 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
+ std::vector< std::vector >
+ phasePressures(const UnstructuredGrid& G,
+ const Region& reg,
+ const CellRange& cells,
+ const double grav)
+ {
+ std::array span =
+ {{ std::numeric_limits::max() ,
+ -std::numeric_limits::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 pval;
+ std::vector 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
+ std::vector< std::vector >
+ phaseSaturations(const Region& reg,
+ const CellRange& cells,
+ const BlackoilPropertiesInterface& props,
+ std::vector< std::vector >& 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 > 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::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
+ std::vector computeRs(const UnstructuredGrid& grid,
+ const CellRangeType& cells,
+ const std::vector oil_pressure,
+ const Miscibility::RsFunction& rs_func,
+ const std::vector gas_saturation)
+ {
+ assert(grid.dimensions == 3);
+ std::vector 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 convertSats(const std::vector< std::vector >& sat)
+ {
+ const int np = sat.size();
+ const int nc = sat[0].size();
+ std::vector 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 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 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
diff --git a/tests/capillary.DATA b/tests/capillary.DATA
new file mode 100644
index 000000000..89c938999
--- /dev/null
+++ b/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
+/
diff --git a/tests/capillary_overlap.DATA b/tests/capillary_overlap.DATA
new file mode 100644
index 000000000..b59e42202
--- /dev/null
+++ b/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
diff --git a/tests/deadfluids.DATA b/tests/deadfluids.DATA
new file mode 100644
index 000000000..89232054b
--- /dev/null
+++ b/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
+/
diff --git a/tests/equil_livegas.DATA b/tests/equil_livegas.DATA
new file mode 100644
index 000000000..a0676f52c
--- /dev/null
+++ b/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
diff --git a/tests/equil_liveoil.DATA b/tests/equil_liveoil.DATA
new file mode 100644
index 000000000..e1b417c46
--- /dev/null
+++ b/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
diff --git a/tests/equil_liveoil_grid.DATA b/tests/equil_liveoil_grid.DATA
new file mode 100644
index 000000000..92e503c5c
--- /dev/null
+++ b/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
+/
diff --git a/tests/equil_rsvd_and_rvvd.DATA b/tests/equil_rsvd_and_rvvd.DATA
new file mode 100644
index 000000000..3076524e3
--- /dev/null
+++ b/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
diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp
new file mode 100644
index 000000000..d68470b08
--- /dev/null
+++ b/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
+#include
+
+/* --- our own headers --- */
+
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+BOOST_AUTO_TEST_SUITE ()
+
+BOOST_AUTO_TEST_CASE (PhasePressure)
+{
+ typedef std::vector PVal;
+ typedef std::vector PPress;
+
+ std::shared_ptr
+ 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 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
+ region(record, calc,
+ std::make_shared(),
+ std::make_shared(),
+ props.phaseUsage());
+
+ std::vector 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 PVal;
+ typedef std::vector PPress;
+
+ std::shared_ptr
+ 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 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 region[] =
+ {
+ Opm::Equil::EquilReg(record[0], calc,
+ std::make_shared(),
+ std::make_shared(),
+ props.phaseUsage())
+ ,
+ Opm::Equil::EquilReg(record[0], calc,
+ std::make_shared(),
+ std::make_shared(),
+ props.phaseUsage())
+ ,
+ Opm::Equil::EquilReg(record[1], calc,
+ std::make_shared(),
+ std::make_shared(),
+ props.phaseUsage())
+ ,
+ Opm::Equil::EquilReg(record[1], calc,
+ std::make_shared(),
+ std::make_shared(),
+ 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 > 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 >::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::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