diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp
new file mode 100644
index 000000000..61f7debc9
--- /dev/null
+++ b/opm/core/simulator/initStateEquil.hpp
@@ -0,0 +1,173 @@
+/*
+ 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
+
+struct UnstructuredGrid;
+
+namespace Opm
+{
+ namespace equil {
+ template
+ class DensityCalculator;
+
+ template <>
+ class DensityCalculator< BlackoilPropertiesInterface > {
+ public:
+ DensityCalculator(const BlackoilPropertiesInterface& props,
+ const int c)
+ : props_(props)
+ , c_(1, c)
+ {
+ }
+
+ 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_;
+ };
+
+ namespace miscibility {
+ struct NoMixing {
+ double
+ operator()(const double /* depth */,
+ const double /* press */) const
+ {
+ return 0.0;
+ }
+ };
+
+ class RsVD {
+ public:
+ RsVD(const std::vector& depth,
+ const std::vector& rs)
+ : depth_(depth)
+ , rs_(rs)
+ {
+ }
+
+ double
+ operator()(const double depth,
+ const double /* press */) const
+ {
+ return linearInterpolation(depth_, rs_, depth);
+ }
+
+ private:
+ std::vector depth_;
+ std::vector rs_;
+ };
+ } // namespace miscibility
+
+ struct EquilRecord {
+ struct {
+ double depth;
+ double press;
+ } main, woc, goc;
+ };
+
+ template
+ class EquilReg {
+ public:
+ EquilReg(const EquilRecord& rec,
+ const DensCalc& density,
+ const RS& rs,
+ const RV& rv,
+ const PhaseUsage& pu)
+ : rec_ (rec)
+ , density_(density)
+ , rs_ (rs)
+ , rv_ (rv)
+ , pu_ (pu)
+ {
+ }
+
+ typedef DensCalc CalcDensity;
+ typedef RS CalcDissolution;
+ typedef RV CalcEvaporation;
+
+ double datum() const { return this->rec_.main.depth; }
+ double pressure() const { return this->rec_.main.press; }
+
+ double zwoc() const { return this->rec_.woc .depth; }
+ double pcow_woc() const { return this->rec_.woc .press; }
+
+ double zgoc() const { return this->rec_.goc .depth; }
+ double pcgo_goc() const { return this->rec_.goc .press; }
+
+ const CalcDensity&
+ densityCalculator() const { return this->density_; }
+
+ const CalcDissolution&
+ dissolutionCalculator() const { return this->rs_; }
+
+ const CalcEvaporation&
+ evaporationCalculator() const { return this->rv_; }
+
+ const PhaseUsage&
+ phaseUsage() const { return this->pu_; }
+
+ private:
+ EquilRecord rec_;
+ DensCalc density_;
+ RS rs_;
+ RV rv_;
+ PhaseUsage pu_;
+ };
+
+ template
+ std::vector< std::vector >
+ phasePressures(const UnstructuredGrid& G,
+ const Region& reg,
+ const double grav = unit::gravity);
+ } // 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..cb19d0205
--- /dev/null
+++ b/opm/core/simulator/initStateEquil_impl.hpp
@@ -0,0 +1,513 @@
+/*
+ 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 std::binary_function {
+ 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();
+ const int i = (x - span_[0]) / h;
+ const double t = (x - (span_[0] + i*h)) / h;
+
+ 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,
+ std::vector& p )
+ {
+ const int nd = G.dimensions, nc = G.number_of_cells;
+ const double* depth = & G.cell_centroids[0*nd + (nd - 1)];
+
+ enum { up = 0, down = 1 };
+
+ for (int c = 0; c < nc; ++c, depth += nd) {
+ const double z = *depth;
+ 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,
+ 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, press);
+ }
+
+ template
+ void
+ oil(const UnstructuredGrid& G ,
+ const Region& reg ,
+ const std::array& span ,
+ const double grav ,
+ 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, 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,
+ 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, press);
+ }
+ } // namespace PhasePressure
+
+ template
+ void
+ equilibrateOWG(const UnstructuredGrid& G,
+ const Region& reg,
+ const double grav,
+ const std::array& span,
+ 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,
+ press[ oix ], po_woc, po_goc);
+ }
+
+ if (PhaseUsed::water(pu)) {
+ const int wix = PhaseIndex::water(pu);
+ PhasePressure::water(G, reg, span, grav,
+ po_woc, press[ wix ]);
+ }
+
+ if (PhaseUsed::gas(pu)) {
+ const int gix = PhaseIndex::gas(pu);
+ PhasePressure::gas(G, reg, span, grav,
+ po_goc, press[ gix ]);
+ }
+ }
+ } // namespace Details
+
+ namespace equil {
+ template
+ std::vector< std::vector >
+ phasePressures(const UnstructuredGrid& G,
+ const Region& reg,
+ const double grav)
+ {
+ std::array span =
+ {{ std::numeric_limits::max() ,
+ -std::numeric_limits::max() }}; // Symm. about 0.
+ {
+ // This code is only supported in three space dimensions
+ assert (G.dimensions == 3);
+
+ const int nd = G.dimensions;
+ const double* depth = & G.node_coordinates[0*nd + (nd - 1)];
+
+ for (int n = 0; n < G.number_of_nodes; ++n, depth += nd) {
+ const double z = *depth;
+
+ 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(G.number_of_cells, 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, press);
+ }
+
+ return press;
+ }
+ } // namespace equil
+} // namespace Opm
+
+#endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp
new file mode 100644
index 000000000..4d7cbe836
--- /dev/null
+++ b/tests/test_equil.cpp
@@ -0,0 +1,89 @@
+/*
+ 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
+
+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,
+ Opm::equil::miscibility::NoMixing(),
+ Opm::equil::miscibility::NoMixing(),
+ props.phaseUsage());
+
+ const double grav = 10;
+ const PPress ppress = Opm::equil::phasePressures(*G, region, 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_SUITE_END()