From 4c39a8a5956442932b96bc452412a11c3fd556bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 14 Jan 2014 20:37:28 +0100 Subject: [PATCH 01/53] Add basic equilibration facility This commit adds a simple facility for calculating initial phase pressures assuming stationary conditions, a known reference pressure in the oil zone as well as the depth and capillary pressures at the water-oil and gas-oil contacts. Function 'Opm::equil::phasePressures()' uses a simple ODE/IVP-based approach, solved using the traditional RK4 method with constant step sizes, to derive the required pressure values. Specifically, we solve the ODE dp/dz = rho(z,p) * g with 'z' represening depth, 'p' being a phase pressure and 'rho' the associate phase density. Finally, 'g' is the acceleration of gravity. We assume that we can calculate phase densities, e.g., from table look-up. This assumption holds in the case of an ECLIPSE input deck. Using RK4 with constant step sizes is a limitation of this implementation. This, basically, assumes that the phase densities varies only smoothly with depth and pressure (at reservoir conditions). --- opm/core/simulator/initStateEquil.hpp | 173 +++++++ opm/core/simulator/initStateEquil_impl.hpp | 513 +++++++++++++++++++++ tests/test_equil.cpp | 89 ++++ 3 files changed, 775 insertions(+) create mode 100644 opm/core/simulator/initStateEquil.hpp create mode 100644 opm/core/simulator/initStateEquil_impl.hpp create mode 100644 tests/test_equil.cpp 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() From c582d00fb6f97ad92daf71b898dd90e541ae6988 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 17 Jan 2014 17:43:27 +0100 Subject: [PATCH 02/53] Compute phase pressures in subset of cells This commit adds support for assigning the initial phase pressure distribution to a subset of the total grid cells. This is needed in order to fully support equilibration regions. The existing region support (template parameter 'Region' in function 'phasePressures()') was only used/needed to define PVT property (specifically, the fluid phase density) calculator pertaining to a particular equilibration region. --- opm/core/simulator/initStateEquil.hpp | 3 +- opm/core/simulator/initStateEquil_impl.hpp | 107 ++++++++++++++++----- tests/test_equil.cpp | 5 +- 3 files changed, 87 insertions(+), 28 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 61f7debc9..b81183cfe 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -160,10 +160,11 @@ namespace Opm PhaseUsage pu_; }; - template + template std::vector< std::vector > phasePressures(const UnstructuredGrid& G, const Region& reg, + const CellRange& cells, const double grav = unit::gravity); } // namespace equil } // namespace Opm diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index cb19d0205..befbf2757 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -286,31 +286,40 @@ namespace Opm } // namespace PhaseIndex namespace PhasePressure { - template + template void assign(const UnstructuredGrid& G , const std::array& f , const double split, + const CellRange& cells, std::vector& p ) { - const int nd = G.dimensions, nc = G.number_of_cells; - const double* depth = & G.cell_centroids[0*nd + (nd - 1)]; + const int nd = G.dimensions; enum { up = 0, down = 1 }; - for (int c = 0; c < nc; ++c, depth += nd) { - const double z = *depth; + 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 + 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; @@ -336,15 +345,17 @@ namespace Opm } }; - assign(G, wpress, z0, press); + assign(G, wpress, z0, cells, press); } - template + 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) @@ -376,7 +387,7 @@ namespace Opm } }; - assign(G, opress, z0, press); + assign(G, opress, z0, cells, press); po_woc = -std::numeric_limits::max(); po_goc = -std::numeric_limits::max(); @@ -394,13 +405,15 @@ namespace Opm else { po_goc = p0; } // GOC *at* datum } - template + 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; @@ -431,16 +444,18 @@ namespace Opm } }; - assign(G, gpress, z0, press); + assign(G, gpress, z0, cells, press); } } // namespace PhasePressure - template + 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(); @@ -448,53 +463,93 @@ namespace Opm double po_woc = -1, po_goc = -1; if (PhaseUsed::oil(pu)) { const int oix = PhaseIndex::oil(pu); - PhasePressure::oil(G, reg, span, grav, + 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, press[ wix ]); + 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, press[ gix ]); + PhasePressure::gas(G, reg, span, grav, po_goc, + cells, press[ gix ]); } } } // namespace Details namespace equil { - template + 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; - const double* depth = & G.node_coordinates[0*nd + (nd - 1)]; + const int nd = G.dimensions; - for (int n = 0; n < G.number_of_nodes; ++n, depth += nd) { - const double z = *depth; + // Define short-name aliases to reduce visual clutter. + const double* const nc = & G.node_coordinates[0]; - if (z < span[0]) { span[0] = z; } - if (z > span[1]) { span[1] = z; } + 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(G.number_of_cells, 0.0)); + std::vector press(np, pval(ncell, 0.0)); const double z0 = reg.datum(); const double zwoc = reg.zwoc (); @@ -502,7 +557,7 @@ namespace Opm if (! ((zgoc > z0) || (z0 > zwoc))) { // Datum depth in oil zone (zgoc <= z0 <= zwoc) - Details::equilibrateOWG(G, reg, grav, span, press); + Details::equilibrateOWG(G, reg, grav, span, cells, press); } return press; diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 4d7cbe836..5a97c74ce 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -75,8 +75,11 @@ BOOST_AUTO_TEST_CASE (PhasePressure) Opm::equil::miscibility::NoMixing(), 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, grav); + 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; From 4114da26bd4af4a54c7a22db9366f4f671b94f2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 17 Jan 2014 19:41:22 +0100 Subject: [PATCH 03/53] Test cell subset phase pressure assignment. --- tests/test_equil.cpp | 112 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 112 insertions(+) diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 5a97c74ce..36d1d4527 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -89,4 +89,116 @@ BOOST_AUTO_TEST_CASE (PhasePressure) 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, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + , + Opm::equil::EquilReg(record[0], calc, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + , + Opm::equil::EquilReg(record[1], calc, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + , + Opm::equil::EquilReg(record[1], calc, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + }; + + const int cdim[] = { 2, 1, 2 }; + int ncoarse = cdim[0]; + for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; } + + std::vector< std::vector > 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_SUITE_END() From 6814b97698392f11bd10ebc83118634c1caa3467 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 17 Jan 2014 20:07:51 +0100 Subject: [PATCH 04/53] Add reverse look-up mapping for region vectors Class RegionMapping<> provides an easy way of extracting the cells that belong to any identified region (e.g., as defined by EQLNUM) of the deck. --- opm/core/simulator/initStateEquil.hpp | 91 +++++++++++++++++++++ tests/test_equil.cpp | 110 ++++++++++++++++++++++++++ 2 files changed, 201 insertions(+) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index b81183cfe..5df47a0c2 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -27,6 +27,7 @@ #include #include +#include #include struct UnstructuredGrid; @@ -102,6 +103,96 @@ namespace Opm }; } // namespace miscibility + template < class Region = std::vector > + class RegionMapping { + public: + explicit + RegionMapping(const Region& reg) + : reg_(reg) + { + rev_.init(reg_); + } + + typedef typename Region::value_type RegionId; + typedef typename Region::size_type CellId; + typedef typename std::vector::const_iterator CellIter; + + class CellRange { + public: + CellRange(const CellIter b, + const CellIter e) + : b_(b), e_(e) + {} + + typedef CellIter iterator; + typedef CellIter const_iterator; + + iterator begin() const { return b_; } + iterator end() const { return e_; } + + private: + iterator b_; + iterator e_; + }; + + RegionId + numRegions() const { return RegionId(rev_.p.size()) - 1; } + + RegionId + region(const CellId c) const { return reg_[c]; } + + CellRange + cells(const RegionId r) const { + const RegionId i = r - rev_.low; + return CellRange(rev_.c.begin() + rev_.p[i + 0], + rev_.c.begin() + rev_.p[i + 1]); + } + + private: + Region reg_; + + struct { + typedef typename std::vector::size_type Pos; + std::vector p; + std::vector c; + RegionId low; + + void + init(const Region& reg) + { + typedef typename Region::const_iterator CI; + const std::pair + m = std::minmax_element(reg.begin(), reg.end()); + + low = *m.first; + + const typename Region::size_type + n = *m.second - low + 1; + + p.resize(n + 1); std::fill(p.begin(), p.end(), Pos(0)); + for (CellId i = 0, nc = reg.size(); i < nc; ++i) { + p[ reg[i] - low + 1 ] += 1; + } + + for (typename std::vector::size_type + i = 1, sz = p.size(); i < sz; ++i) { + p[0] += p[i]; + p[i] = p[0] - p[i]; + } + + assert (p[0] == + static_cast(reg.size())); + + c.resize(reg.size()); + for (CellId i = 0, nc = reg.size(); i < nc; ++i) { + c[ p[ reg[i] - low + 1 ] ++ ] = i; + } + + p[0] = 0; + } + } rev_; + }; + struct EquilRecord { struct { double depth; diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 36d1d4527..80ae26529 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -25,6 +25,8 @@ #include #include +#include + #include #include @@ -201,4 +203,112 @@ BOOST_AUTO_TEST_CASE (CellSubset) BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol); } +BOOST_AUTO_TEST_CASE (RegMapping) +{ + 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, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + , + Opm::equil::EquilReg(record[0], calc, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + , + Opm::equil::EquilReg(record[1], calc, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + , + Opm::equil::EquilReg(record[1], calc, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + }; + + std::vector eqlnum(G->number_of_cells); + { + std::vector cells(G->number_of_cells); + std::iota(cells.begin(), cells.end(), 0); + + const int cdim[] = { 2, 1, 2 }; + int ncoarse = cdim[0]; + for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; } + + partition_unif_idx(G->dimensions, G->number_of_cells, + G->cartdims, cdim, + &cells[0], &eqlnum[0]); + } + Opm::equil::RegionMapping<> eqlmap(eqlnum); + + PPress ppress(2, PVal(G->number_of_cells, 0)); + for (int r = 0, e = eqlmap.numRegions(); r != e; ++r) + { + const Opm::equil::RegionMapping<>::CellRange& + rng = eqlmap.cells(r); + + const int rno = r; + const double grav = 10; + const PPress p = + Opm::equil::phasePressures(*G, region[rno], rng, grav); + + PVal::size_type i = 0; + for (Opm::equil::RegionMapping<>::CellRange::const_iterator + c = rng.begin(), ce = rng.end(); + c != ce; ++c, ++i) + { + assert (i < p[0].size()); + + ppress[0][*c] = p[0][i]; + ppress[1][*c] = p[1][i]; + } + } + + const int first = 0, last = G->number_of_cells - 1; + const double reltol = 1.0e-8; + BOOST_CHECK_CLOSE(ppress[0][first] , 105e3 , reltol); + BOOST_CHECK_CLOSE(ppress[0][last ] , 195e3 , reltol); + BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol); + BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol); +} + BOOST_AUTO_TEST_SUITE_END() From 955cb795f3ccc0cab509f67ac2d7a22e103145a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 20 Jan 2014 01:25:33 +0100 Subject: [PATCH 05/53] Document public interface of phasePressures(). --- opm/core/simulator/initStateEquil.hpp | 356 ++++++++++++++++++++++++-- 1 file changed, 336 insertions(+), 20 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 5df47a0c2..43f77aac4 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -30,17 +30,45 @@ #include #include +/** + * \file + * Facilities for an ECLIPSE-style equilibration-based + * initialisation scheme (keyword 'EQUIL'). + */ struct UnstructuredGrid; 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) @@ -48,6 +76,16 @@ namespace Opm { } + /** + * 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 @@ -71,8 +109,28 @@ namespace Opm const std::vector c_; }; + /** + * Types and routines relating to phase mixing in + * equilibration calculations. + */ namespace miscibility { + /** + * Type that implements "no phase mixing" policy. + */ struct NoMixing { + /** + * 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 @@ -81,8 +139,19 @@ namespace Opm } }; + /** + * Type that implements "dissolved gas-oil ratio" + * tabulated as a function of depth policy. Data + * typically taken from keyword 'RSVD'. + */ class RsVD { public: + /** + * Constructor. + * + * \param[in] depth Depth nodes. + * \param[in] rs Dissolved gas-oil ratio at @c depth. + */ RsVD(const std::vector& depth, const std::vector& rs) : depth_(depth) @@ -90,6 +159,18 @@ namespace Opm { } + /** + * 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 @@ -98,14 +179,31 @@ namespace Opm } private: - std::vector depth_; - std::vector rs_; + std::vector depth_; /**< Depth nodes */ + std::vector rs_; /**< Dissolved gas-oil ratio */ }; } // namespace miscibility + /** + * Forward and reverse mappings between cells and + * regions/partitions (e.g., the ECLIPSE-style 'SATNUM', + * 'PVTNUM', or 'EQUILNUM' arrays). + * + * \tparam Region Type of a forward region mapping. Expected + * to provide indexed access through + * operator[]() as well as inner types + * 'value_type', 'size_type', and + * 'const_iterator'. + */ template < class Region = std::vector > class RegionMapping { public: + /** + * Constructor. + * + * \param[in] reg Forward region mapping, restricted to + * active cells only. + */ explicit RegionMapping(const Region& reg) : reg_(reg) @@ -113,34 +211,82 @@ namespace Opm rev_.init(reg_); } - typedef typename Region::value_type RegionId; - typedef typename Region::size_type CellId; + /** + * Type of forward (cell-to-region) mapping result. + * Expected to be an integer. + */ + typedef typename Region::value_type RegionId; + + /** + * Type of reverse (region-to-cell) mapping (element) + * result. + */ + typedef typename Region::size_type CellId; + + /** + * Type of reverse region-to-cell range bounds and + * iterators. + */ typedef typename std::vector::const_iterator CellIter; + /** + * Range of cells. Result from reverse (region-to-cell) + * mapping. + */ class CellRange { public: + /** + * Constructor. + * + * \param[in] b Beginning of range. + * \param[in] e One past end of range. + */ CellRange(const CellIter b, const CellIter e) : b_(b), e_(e) {} - typedef CellIter iterator; + /** + * Read-only iterator on cell ranges. + */ typedef CellIter const_iterator; - iterator begin() const { return b_; } - iterator end() const { return e_; } + /** + * Beginning of cell range. + */ + const_iterator begin() const { return b_; } + + /** + * One past end of cell range. + */ + const_iterator end() const { return e_; } private: - iterator b_; - iterator e_; + const_iterator b_; + const_iterator e_; }; + /** + * Number of declared regions in cell-to-region mapping. + */ RegionId numRegions() const { return RegionId(rev_.p.size()) - 1; } + /** + * Compute region number of given active cell. + * + * \param[in] c Active cell + * \return Region to which @c c belongs. + */ RegionId region(const CellId c) const { return reg_[c]; } + /** + * Extract active cells in particular region. + * + * \param[in] r Region number + * \returns Range of active cells in region @c r. + */ CellRange cells(const RegionId r) const { const RegionId i = r - rev_.low; @@ -149,14 +295,24 @@ namespace Opm } private: + /** + * Copy of forward region mapping (cell-to-region). + */ Region reg_; + /** + * Reverse mapping (region-to-cell). + */ struct { typedef typename std::vector::size_type Pos; - std::vector p; - std::vector c; - RegionId low; + std::vector p; /**< Region start pointers */ + std::vector c; /**< Region cells */ + RegionId low; /**< Smallest region number */ + /** + * Compute reverse mapping. Standard linear insertion + * sort algorithm. + */ void init(const Region& reg) { @@ -190,9 +346,32 @@ namespace Opm p[0] = 0; } - } rev_; + } rev_; /**< Reverse mapping instance */ }; + /** + * 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; @@ -200,11 +379,63 @@ namespace Opm } main, woc, goc; }; + /** + * 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. + * + * \tparam RS Type that provides access to a calculator for + * (initial) dissolved gas-oil ratios as a function of depth + * and (oil) pressure. Must implement an operator() declared + * as + * + * double + * operator()(const double depth, + * const double press) + * + * that calculates the dissolved gas-oil ratio at depth @c + * depth and (oil) pressure @c press. + * + * \tparam RV Type that provides access to a calculator for + * (initial) vapourised oil-gas ratios as a function of depth + * and (gas) pressure. Must implement an operator() declared + * as + * + * double + * operator()(const double depth, + * const double press) + * + * that calculates the vapourised oil-gas ratio at depth @c + * depth and (gas) 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, const RS& rs, @@ -218,39 +449,124 @@ namespace Opm { } + /** + * Type of density calculator. + */ typedef DensCalc CalcDensity; - typedef RS CalcDissolution; - typedef RV CalcEvaporation; + /** + * Type of dissolved gas-oil ratio calculator. + */ + typedef RS CalcDissolution; + + /** + * Type of vapourised oil-gas ratio calculator. + */ + typedef RV 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_; - DensCalc density_; - RS rs_; - RV rv_; - PhaseUsage pu_; + EquilRecord rec_; /**< Equilibration data */ + DensCalc density_; /**< Density calculator */ + RS rs_; /**< RS calculator */ + RV rv_; /**< RV calculator */ + PhaseUsage pu_; /**< Active phase summary */ }; + /** + * 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. + * + * \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, From ada48fd0b5b363fd14152a345cfecf789638d039 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 20 Jan 2014 10:33:42 +0100 Subject: [PATCH 06/53] Document requirements of CellRange. --- opm/core/simulator/initStateEquil.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 43f77aac4..4d7a6f0ff 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -555,7 +555,10 @@ namespace Opm * * \tparam CellRange Type of cell range that demarcates the * cells pertaining to the current - * equilibration region. + * 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. From 76844c9c59ba720d4de2bf428f33e5c56807ab81 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 21 Jan 2014 13:49:06 +0100 Subject: [PATCH 07/53] Include for std::iota() Header was missing in earlier revision. --- tests/test_equil.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 80ae26529..27b46b559 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include From f71ba5a9f7fd014194ae6b1f8769a1bbe3575c29 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 21 Jan 2014 17:49:02 +0100 Subject: [PATCH 08/53] Install crude handling of data point outside vertical span The initial implementation of RK4IVP<>::operator() failed to take into account the possibility that we might need to evaluate the function outside the vertical span for which it was initially defined. This situation occurs, for instance, in the not uncommon cases of the GOC being above or the WOC being below the model. This commit installs a crude Hermitian extrapolation procedure to handle these cases. Refinements are likely. --- opm/core/simulator/initStateEquil_impl.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index befbf2757..ffc566be3 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -73,9 +73,13 @@ namespace Opm // Dense output (O(h**3)) according to Shampine // (Hermite interpolation) const double h = stepsize(); - const int i = (x - span_[0]) / h; + 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]; From 311ac78340a8895535a59d3de457b758e1d16fa1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 23 Jan 2014 10:16:49 +0100 Subject: [PATCH 09/53] Add a layer of glue to extract data from deck This is a work in progress. --- opm/core/simulator/initStateEquil.hpp | 126 ++++++++++++++++++++++++++ 1 file changed, 126 insertions(+) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 4d7a6f0ff..f9b5a65af 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -20,6 +20,7 @@ #ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED #define OPM_INITSTATEEQUIL_HEADER_INCLUDED +#include #include #include #include @@ -576,6 +577,131 @@ namespace Opm const Region& reg, const CellRange& cells, const double grav = unit::gravity); + + 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_ } + }; + + 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")) { + eqlnum = deck.getIntegerValue("EQLNUM"); + } + else { + // No explicit equilibration region. + // All cells in region zero. + eqlnum.assign(G.number_of_cells, 0); + } + + return eqlnum; + } + + template + class PhasePressureComputer; + + template <> + class PhasePressureComputer { + public: + PhasePressureComputer(const BlackoilPropertiesInterface& props, + const EclipseGridParser& deck , + const UnstructuredGrid& G ) + : pp_(props.numPhases(), + std::vector(G.number_of_cells)) + { + const std::vector rec = getEquil(deck); + const RegionMapping<> eqlmap(equilnum(deck, G)); + + calcII(eqlmap, rec, props, G); + } + + typedef std::vector PVal; + typedef std::vector PPress; + + const PPress& press() const { return pp_; } + + private: + typedef DensityCalculator RhoCalc; + typedef EquilReg EqReg; + + PPress pp_; + + template + void + calcII(const RMap& reg , + const std::vector< EquilRecord >& rec , + const Opm::BlackoilPropertiesInterface& props, + const UnstructuredGrid& G ) + { + 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, NoMix(), NoMix(), + props.phaseUsage()); + + const PPress& res = phasePressures(G, eqreg, cells); + + for (int p = 0, np = props.numPhases(); p < np; ++p) { + PVal& d = pp_[p]; + PVal::const_iterator s = res[p].begin(); + for (typename RMap::CellRange::const_iterator + c = cells.begin(), + e = cells.end(); + c != e; ++c, ++s) + { + d[*c] = *s; + } + } + } + } + }; + } // namespace DeckDependent } // namespace equil } // namespace Opm From 30d65e6c32a267b689a8af589d13f4e9d50700ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 3 Feb 2014 11:32:46 +0100 Subject: [PATCH 10/53] Removed RK4IVP's inheritance from binary_function. Three reasons: - class is a unary functor, - the typedefs obtained were not used, - binary_function is deprecated in C++11. --- opm/core/simulator/initStateEquil_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index ffc566be3..381a3f51c 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -32,7 +32,7 @@ namespace Opm { namespace Details { template - class RK4IVP : public std::binary_function { + class RK4IVP { public: RK4IVP(const RHS& f , const std::array& span, From 85be4b1263cb0899acc89b500ec43adacb19556f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 3 Feb 2014 15:36:20 +0100 Subject: [PATCH 11/53] Add (unfinished) test case. --- tests/deadfluids.DATA | 9 +++++++++ tests/test_equil.cpp | 20 ++++++++++++++++++++ 2 files changed, 29 insertions(+) create mode 100644 tests/deadfluids.DATA diff --git a/tests/deadfluids.DATA b/tests/deadfluids.DATA new file mode 100644 index 000000000..37b3c6ee4 --- /dev/null +++ b/tests/deadfluids.DATA @@ -0,0 +1,9 @@ +WATER +OIL +GAS + +PVDO +100 1.0 1.0 +200 0.9 1.0 +/ + diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 27b46b559..7060c501c 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -23,6 +23,7 @@ #include #include +#include #include #include @@ -92,6 +93,9 @@ BOOST_AUTO_TEST_CASE (PhasePressure) BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol); } + + + BOOST_AUTO_TEST_CASE (CellSubset) { typedef std::vector PVal; @@ -204,6 +208,9 @@ BOOST_AUTO_TEST_CASE (CellSubset) BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol); } + + + BOOST_AUTO_TEST_CASE (RegMapping) { typedef std::vector PVal; @@ -312,4 +319,17 @@ BOOST_AUTO_TEST_CASE (RegMapping) BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol); } + + +BOOST_AUTO_TEST_CASE (DeckAllDead) +{ + std::shared_ptr + grid(create_grid_cart3d(10, 1, 10), destroy_grid); + Opm::EclipseGridParser deck("deadfluids.DATA"); + Opm::BlackoilPropertiesFromDeck props(deck, *grid, false); + Opm::equil::DeckDependent::PhasePressureComputer comp(props, deck, *grid); +} + + + BOOST_AUTO_TEST_SUITE_END() From 634b78a9d99ab8d12eff2a944bd4e8da85204ae6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 4 Feb 2014 09:38:47 +0100 Subject: [PATCH 12/53] Created simple data for init testing. --- tests/deadfluids.DATA | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/tests/deadfluids.DATA b/tests/deadfluids.DATA index 37b3c6ee4..190d5d975 100644 --- a/tests/deadfluids.DATA +++ b/tests/deadfluids.DATA @@ -7,3 +7,25 @@ PVDO 200 0.9 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 2 0 8 0 +/ From 45de38a01982011e17752e6af921c27421879b5f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 4 Feb 2014 10:56:09 +0100 Subject: [PATCH 13/53] Throw exception if datum not in oil zone. We are not capable of handling this, and must abort. --- opm/core/simulator/initStateEquil_impl.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index 381a3f51c..a713e3b34 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -562,6 +562,8 @@ namespace Opm 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; From 04eea929bac27f4879e28a6a7007917b1cc98741 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 4 Feb 2014 10:57:23 +0100 Subject: [PATCH 14/53] Fix contact depths in test deck. --- tests/deadfluids.DATA | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/deadfluids.DATA b/tests/deadfluids.DATA index 190d5d975..8c368ff88 100644 --- a/tests/deadfluids.DATA +++ b/tests/deadfluids.DATA @@ -27,5 +27,5 @@ DENSITY / EQUIL -5 150 2 0 8 0 +5 150 8 0 2 0 / From b134fe048ad10c5c28a945b18e8e5cd36a75b490 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 5 Feb 2014 11:26:29 +0100 Subject: [PATCH 15/53] Still working on test_equil.cpp. --- tests/test_equil.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 7060c501c..844cf6ad2 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -324,10 +324,19 @@ BOOST_AUTO_TEST_CASE (RegMapping) BOOST_AUTO_TEST_CASE (DeckAllDead) { std::shared_ptr - grid(create_grid_cart3d(10, 1, 10), destroy_grid); + grid(create_grid_cart3d(1, 1, 10), destroy_grid); Opm::EclipseGridParser deck("deadfluids.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, *grid, false); Opm::equil::DeckDependent::PhasePressureComputer comp(props, deck, *grid); + const auto& pressures = comp.press(); + BOOST_REQUIRE(pressures.size() == 3); + BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells); + for (auto pp : pressures) { + for (auto p : pp){ + std::cout << p << ' '; + } + std::cout << std::endl; + } } From 82216fa24f628a617f43723d60cbc41fc4e2eee1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 19 Feb 2014 13:38:21 +0100 Subject: [PATCH 16/53] Add (defaulted) gravity argument in some places. This is done to facilitate testing, using gravity = 10 m/s^2 for example. --- opm/core/simulator/initStateEquil.hpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index f9b5a65af..cde37352b 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -643,14 +643,15 @@ namespace Opm public: PhasePressureComputer(const BlackoilPropertiesInterface& props, const EclipseGridParser& deck , - const UnstructuredGrid& G ) + const UnstructuredGrid& G , + const double grav = unit::gravity) : pp_(props.numPhases(), std::vector(G.number_of_cells)) { const std::vector rec = getEquil(deck); const RegionMapping<> eqlmap(equilnum(deck, G)); - calcII(eqlmap, rec, props, G); + calcII(eqlmap, rec, props, G, grav); } typedef std::vector PVal; @@ -669,7 +670,8 @@ namespace Opm calcII(const RMap& reg , const std::vector< EquilRecord >& rec , const Opm::BlackoilPropertiesInterface& props, - const UnstructuredGrid& G ) + const UnstructuredGrid& G , + const double grav) { typedef miscibility::NoMixing NoMix; @@ -685,7 +687,7 @@ namespace Opm const EqReg eqreg(rec[r], calc, NoMix(), NoMix(), props.phaseUsage()); - const PPress& res = phasePressures(G, eqreg, cells); + const PPress& res = phasePressures(G, eqreg, cells, grav); for (int p = 0, np = props.numPhases(); p < np; ++p) { PVal& d = pp_[p]; From 2ebaef62a00cfbc52f36cea4e1b5144718571671 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 19 Feb 2014 13:40:02 +0100 Subject: [PATCH 17/53] Modify test data. --- tests/deadfluids.DATA | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/deadfluids.DATA b/tests/deadfluids.DATA index 8c368ff88..ffaebc7ff 100644 --- a/tests/deadfluids.DATA +++ b/tests/deadfluids.DATA @@ -4,7 +4,7 @@ GAS PVDO 100 1.0 1.0 -200 0.9 1.0 +200 0.5 1.0 / PVDG @@ -27,5 +27,5 @@ DENSITY / EQUIL -5 150 8 0 2 0 +5 150 5 0 2 0 / From ebb6eaf74211bb931e552195e5e8d163f29b5294 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 19 Feb 2014 13:41:20 +0100 Subject: [PATCH 18/53] Complete pressure test for dead-oil deck. --- tests/test_equil.cpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 844cf6ad2..ff48cbc1a 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -327,16 +327,20 @@ BOOST_AUTO_TEST_CASE (DeckAllDead) grid(create_grid_cart3d(1, 1, 10), destroy_grid); Opm::EclipseGridParser deck("deadfluids.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, *grid, false); - Opm::equil::DeckDependent::PhasePressureComputer comp(props, deck, *grid); + Opm::equil::DeckDependent::PhasePressureComputer comp(props, deck, *grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells); - for (auto pp : pressures) { - for (auto p : pp){ - std::cout << p << ' '; - } - std::cout << std::endl; - } + + const int first = 0, last = grid->number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-3; + BOOST_CHECK_CLOSE(pressures[0][first] , 14955e3 , reltol); + BOOST_CHECK_CLOSE(pressures[0][last ] , 15045e3 , reltol); + BOOST_CHECK_CLOSE(pressures[1][last] , 1.50473e7 , reltol); } From b09dd7750a1961ae2b123c167de61cca332007ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 19 Feb 2014 13:42:07 +0100 Subject: [PATCH 19/53] Add saturation init facilities. This adds the function phaseSaturations() and some helpers: satFromPc() and satFromSumOfPcs(). --- opm/core/simulator/initStateEquil.hpp | 223 ++++++++++++++++++++++++++ 1 file changed, 223 insertions(+) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index cde37352b..6613d550a 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -533,6 +534,143 @@ namespace Opm 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 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; + } + + + /** * Compute initial phase pressures by means of equilibration. * @@ -578,6 +716,91 @@ namespace Opm 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) + { + 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. + + const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid]; + 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; + // Find saturations from pressure differences by + // inverting capillary pressure functions. + double sw = 0.0; + if (reg.phaseUsage().phase_used[BlackoilPhases::Aqua]) { + const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua]; + 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 (reg.phaseUsage().phase_used[BlackoilPhases::Vapour]) { + const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour]; + // 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; + } + if (sg > 0.0 && sw > 0.0) { + // Overlapping gas-oil and oil-water transition + // zones. Must recalculate using gas-water + // capillary pressure. + const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua]; + const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour]; + 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; + } + phase_saturations[oilpos][local_index] = 1.0 - sw - sg; + } + } + + + namespace DeckDependent { inline std::vector From e1069cf39b6acdf4578cd970bf28582aa49737e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 20 Feb 2014 15:24:27 +0100 Subject: [PATCH 20/53] Fix bugs in saturation initialisation and helpers. --- opm/core/simulator/initStateEquil.hpp | 38 ++++++++++++++++++--------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 6613d550a..723776c14 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -661,12 +661,20 @@ namespace Opm // Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc const PcEqSum f(props, phase1, phase2, cell, target_pc); - 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; + 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; + } } @@ -760,35 +768,38 @@ namespace Opm } 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 (reg.phaseUsage().phase_used[BlackoilPhases::Aqua]) { - const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua]; + 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 (reg.phaseUsage().phase_used[BlackoilPhases::Vapour]) { - const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour]; + 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; } - if (sg > 0.0 && sw > 0.0) { + if (gas && water && sg > smin[gaspos] && sw > smin[waterpos]) { // Overlapping gas-oil and oil-water transition // zones. Must recalculate using gas-water // capillary pressure. - const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua]; - const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour]; const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index]; sw = satFromSumOfPcs(props, waterpos, gaspos, cell, pcgw); sg = 1.0 - sw; @@ -797,6 +808,7 @@ namespace Opm } phase_saturations[oilpos][local_index] = 1.0 - sw - sg; } + return phase_saturations; } From 4ec65f2f8f47da1528e7f576308548d9a9d9057b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 20 Feb 2014 15:39:15 +0100 Subject: [PATCH 21/53] Add another test deck for initialisation. This deck includes capillary functions. --- tests/capillary.DATA | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) create mode 100644 tests/capillary.DATA diff --git a/tests/capillary.DATA b/tests/capillary.DATA new file mode 100644 index 000000000..d8d163fcd --- /dev/null +++ b/tests/capillary.DATA @@ -0,0 +1,31 @@ +WATER +OIL +GAS + +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 +/ From dc6bcead912c8c1804aba29d14d5b44c2ee4da1e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 21 Feb 2014 08:32:15 +0100 Subject: [PATCH 22/53] Add test case for capillary inversion. --- tests/test_equil.cpp | 55 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index ff48cbc1a..712ee31c4 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -21,6 +21,7 @@ #include #include +#include #include #include @@ -345,4 +346,58 @@ BOOST_AUTO_TEST_CASE (DeckAllDead) +BOOST_AUTO_TEST_CASE (CapillaryInversion) +{ + // Test setup. + Opm::GridManager gm(1, 1, 40, 1.0, 1.0, 2.5); + const UnstructuredGrid& grid = *(gm.c_grid()); + Opm::EclipseGridParser deck("capillary.DATA"); + Opm::BlackoilPropertiesFromDeck props(deck, grid, false); + + // Test the capillary inversion for oil-water. + const int cell = 0; + const double reltol = 1.0e-7; + { + const int phase = 0; + const bool increasing = false; + const std::vector pc = { 10.0e5, 0.5e5, 0.4e5, 0.3e5, 0.2e5, 0.1e5, 0.099e5, 0.0e5, -10.0e5 }; + const std::vector s = { 0.2, 0.2, 0.2, 0.466666666666, 0.733333333333, 1.0, 1.0, 1.0, 1.0 }; + BOOST_REQUIRE(pc.size() == s.size()); + for (size_t i = 0; i < pc.size(); ++i) { + const double s_computed = Opm::equil::satFromPc(props, phase, cell, pc[i], increasing); + BOOST_CHECK_CLOSE(s_computed, s[i], reltol); + } + } + + // Test the capillary inversion for gas-oil. + { + const int phase = 2; + const bool increasing = true; + const std::vector pc = { 10.0e5, 0.6e5, 0.5e5, 0.4e5, 0.3e5, 0.2e5, 0.1e5, 0.0e5, -10.0e5 }; + const std::vector s = { 0.8, 0.8, 0.8, 0.533333333333, 0.266666666666, 0.0, 0.0, 0.0, 0.0 }; + BOOST_REQUIRE(pc.size() == s.size()); + for (size_t i = 0; i < pc.size(); ++i) { + const double s_computed = Opm::equil::satFromPc(props, phase, cell, pc[i], increasing); + BOOST_CHECK_CLOSE(s_computed, s[i], reltol); + } + } + + // Test the capillary inversion for gas-water. + { + const int water = 0; + const int gas = 2; + const std::vector pc = { 0.9e5, 0.8e5, 0.6e5, 0.4e5, 0.3e5 }; + const std::vector s = { 0.2, 0.333333333333, 0.6, 0.866666666666, 1.0 }; + BOOST_REQUIRE(pc.size() == s.size()); + for (size_t i = 0; i < pc.size(); ++i) { + const double s_computed = Opm::equil::satFromSumOfPcs(props, water, gas, cell, pc[i]); + BOOST_CHECK_CLOSE(s_computed, s[i], reltol); + } + } +} + + + + + BOOST_AUTO_TEST_SUITE_END() From b2be489e6e048249a6a978c505ed0e07e282e84e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 21 Feb 2014 08:52:25 +0100 Subject: [PATCH 23/53] Add saturation computation to and rename computer class. Opm::equil::DeckDependent::PhasePressureComputer -> Opm::equil::DeckDependent::PhasePressureSaturationComputer --- opm/core/simulator/initStateEquil.hpp | 81 +++++++++++++++++++++++---- tests/test_equil.cpp | 4 +- 2 files changed, 70 insertions(+), 15 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 723776c14..e3460bc97 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -871,42 +871,47 @@ namespace Opm } template - class PhasePressureComputer; + class PhasePressureSaturationComputer; template <> - class PhasePressureComputer { + class PhasePressureSaturationComputer { public: - PhasePressureComputer(const BlackoilPropertiesInterface& props, - const EclipseGridParser& deck , - const UnstructuredGrid& G , - const double grav = unit::gravity) + PhasePressureSaturationComputer(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)) { const std::vector rec = getEquil(deck); const RegionMapping<> eqlmap(equilnum(deck, G)); - calcII(eqlmap, rec, props, G, grav); + calcPressII(eqlmap, rec, props, G, grav); + calcSat(eqlmap, rec, props, G, grav); } typedef std::vector PVal; typedef std::vector PPress; const PPress& press() const { return pp_; } + const PPress& saturation() const { return sat_; } private: typedef DensityCalculator RhoCalc; typedef EquilReg EqReg; PPress pp_; + PPress sat_; template void - calcII(const RMap& reg , - const std::vector< EquilRecord >& rec , - const Opm::BlackoilPropertiesInterface& props, - const UnstructuredGrid& G , - const double grav) + calcPressII(const RMap& reg , + const std::vector< EquilRecord >& rec , + const Opm::BlackoilPropertiesInterface& props, + const UnstructuredGrid& G , + const double grav) { typedef miscibility::NoMixing NoMix; @@ -937,6 +942,58 @@ namespace Opm } } } + + template + void + calcSat(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, NoMix(), NoMix(), + props.phaseUsage()); + + const PPress press = phasePressures(G, eqreg, cells, grav); + const PPress sat = phaseSaturations(eqreg, cells, props, press); + + for (int p = 0, np = props.numPhases(); p < np; ++p) { + PVal& d = pp_[p]; + PVal::const_iterator s = press[p].begin(); + for (typename RMap::CellRange::const_iterator + c = cells.begin(), + e = cells.end(); + c != e; ++c, ++s) + { + d[*c] = *s; + } + } + for (int p = 0, np = props.numPhases(); p < np; ++p) { + PVal& d = sat_[p]; + PVal::const_iterator s = sat[p].begin(); + for (typename RMap::CellRange::const_iterator + c = cells.begin(), + e = cells.end(); + c != e; ++c, ++s) + { + d[*c] = *s; + } + } + } + } + + }; } // namespace DeckDependent } // namespace equil diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 712ee31c4..831482185 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -328,7 +328,7 @@ BOOST_AUTO_TEST_CASE (DeckAllDead) grid(create_grid_cart3d(1, 1, 10), destroy_grid); Opm::EclipseGridParser deck("deadfluids.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, *grid, false); - Opm::equil::DeckDependent::PhasePressureComputer comp(props, deck, *grid, 10.0); + Opm::equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, *grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells); @@ -398,6 +398,4 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion) - - BOOST_AUTO_TEST_SUITE_END() From 7d63cb920471df78eefba6cf7b723d889fa9ff6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 21 Feb 2014 08:55:15 +0100 Subject: [PATCH 24/53] Add test case with capillary transition region. --- tests/test_equil.cpp | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 831482185..d62aa5f1c 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -398,4 +398,42 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion) +BOOST_AUTO_TEST_CASE (DeckWithCapillary) +{ + Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0); + const UnstructuredGrid& grid = *(gm.c_grid()); + Opm::EclipseGridParser deck("capillary.DATA"); + Opm::BlackoilPropertiesFromDeck props(deck, grid, false); + + Opm::equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, grid, 10.0); + const auto& pressures = comp.press(); + BOOST_REQUIRE(pressures.size() == 3); + BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-6; + BOOST_CHECK_CLOSE(pressures[0][first] , 1.45e7 , reltol); + BOOST_CHECK_CLOSE(pressures[0][last ] , 1.545e7 , reltol); + BOOST_CHECK_CLOSE(pressures[1][last] , 1.5351621345e7 , reltol); + + const auto& sats = comp.saturation(); + const std::vector s[3]{ + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.425893333333, 0.774026666666, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0.00736, 0.792746666666, 0.8, 0.8, 0.8, 0.8, 0.574106666666, 0.225973333333, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.79264, 0.007253333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + for (int phase = 0; phase < 3; ++phase) { + BOOST_REQUIRE(sats[phase].size() == s[phase].size()); + for (size_t i = 0; i < s[phase].size(); ++i) { + BOOST_CHECK_CLOSE(sats[phase][i], s[phase][i], reltol); + } + } + +} + + BOOST_AUTO_TEST_SUITE_END() From 6c8babc7d2dddfd5882984aeb3134c8ab91914cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 21 Feb 2014 14:47:14 +0100 Subject: [PATCH 25/53] Fix bug in saturation initialisation. We shall only use gas-water capillary to initialise when we would get unphysical saturations otherwise. --- opm/core/simulator/initStateEquil.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index e3460bc97..8e2b561ff 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -796,9 +796,10 @@ namespace Opm sg = satFromPc(props, gaspos, cell, pcog, increasing); phase_saturations[gaspos][local_index] = sg; } - if (gas && water && sg > smin[gaspos] && sw > smin[waterpos]) { + if (gas && water && (sg + sw > 1.0)) { // Overlapping gas-oil and oil-water transition - // zones. Must recalculate using gas-water + // 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); From 83d04870971864fa1479a6fe5e30da36535f76a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 21 Feb 2014 14:50:45 +0100 Subject: [PATCH 26/53] Add test case with overlapping transitions. Capillary pressure functions and contact depths have been modified to ensure a large overlap. --- tests/capillary_overlap.DATA | 31 +++++++++++++++++++++++++ tests/test_equil.cpp | 44 ++++++++++++++++++++++++++++++++++++ 2 files changed, 75 insertions(+) create mode 100644 tests/capillary_overlap.DATA diff --git a/tests/capillary_overlap.DATA b/tests/capillary_overlap.DATA new file mode 100644 index 000000000..185784eee --- /dev/null +++ b/tests/capillary_overlap.DATA @@ -0,0 +1,31 @@ +WATER +OIL +GAS + +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 +/ + +DENSITY +700 1000 1 +/ + +EQUIL +50 150 50 0.25 45 0.35 +/ diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index d62aa5f1c..c4db1786a 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -432,7 +432,51 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary) BOOST_CHECK_CLOSE(sats[phase][i], s[phase][i], reltol); } } +} + + +BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap) +{ + Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0); + const UnstructuredGrid& grid = *(gm.c_grid()); + Opm::EclipseGridParser deck("capillary_overlap.DATA"); + Opm::BlackoilPropertiesFromDeck props(deck, grid, false); + + Opm::equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, grid, 10.0); + const auto& pressures = comp.press(); + BOOST_REQUIRE(pressures.size() == 3); + BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-6; + BOOST_CHECK_CLOSE(pressures[0][first] , 1.45e7 , reltol); + BOOST_CHECK_CLOSE(pressures[0][last ] , 1.545e7 , reltol); + BOOST_CHECK_CLOSE(pressures[1][last] , 1.5351621345e7 , reltol); + + const auto& sats = comp.saturation(); + // std::cout << "Saturations:\n"; + // for (const auto& sat : sats) { + // for (const double s : sat) { + // std::cout << s << ' '; + // } + // std::cout << std::endl; + // } + const std::vector s[3]{ + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.223141818182, 0.532269090909, 0.78471, 0.91526, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.207743333333, 0.08474, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.776858181818, 0.467730909091, 0.0075466666666, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + for (int phase = 0; phase < 3; ++phase) { + BOOST_REQUIRE(sats[phase].size() == s[phase].size()); + for (size_t i = 0; i < s[phase].size(); ++i) { + BOOST_CHECK_CLOSE(sats[phase][i], s[phase][i], reltol); + } + } } From c689bc757f7dc80179fc6981d5d00cee86160ac3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 24 Feb 2014 13:47:03 +0100 Subject: [PATCH 27/53] Added class RsSatAtContact (not tested). --- opm/core/simulator/initStateEquil.hpp | 74 +++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 8e2b561ff..4ed2b6190 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -184,8 +184,82 @@ namespace Opm std::vector depth_; /**< Depth nodes */ std::vector rs_; /**< Dissolved gas-oil ratio */ }; + + + /** + * 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: + /** + * 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 + { + return std::max(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]; + } + }; } // namespace miscibility + /** * Forward and reverse mappings between cells and * regions/partitions (e.g., the ECLIPSE-style 'SATNUM', From 2b0fcfe7486ff04f8a0e2f97e8c3edec51a7a4a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 24 Feb 2014 15:19:04 +0100 Subject: [PATCH 28/53] Move RegionMapping class to its own header, add test. Class now resides in opm/core/utility/RegionMapping.hpp. --- opm/core/simulator/initStateEquil.hpp | 165 +------------------------- tests/test_equil.cpp | 6 +- 2 files changed, 4 insertions(+), 167 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 4ed2b6190..fd1c6721a 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -260,170 +261,6 @@ namespace Opm } // namespace miscibility - /** - * Forward and reverse mappings between cells and - * regions/partitions (e.g., the ECLIPSE-style 'SATNUM', - * 'PVTNUM', or 'EQUILNUM' arrays). - * - * \tparam Region Type of a forward region mapping. Expected - * to provide indexed access through - * operator[]() as well as inner types - * 'value_type', 'size_type', and - * 'const_iterator'. - */ - template < class Region = std::vector > - class RegionMapping { - public: - /** - * Constructor. - * - * \param[in] reg Forward region mapping, restricted to - * active cells only. - */ - explicit - RegionMapping(const Region& reg) - : reg_(reg) - { - rev_.init(reg_); - } - - /** - * Type of forward (cell-to-region) mapping result. - * Expected to be an integer. - */ - typedef typename Region::value_type RegionId; - - /** - * Type of reverse (region-to-cell) mapping (element) - * result. - */ - typedef typename Region::size_type CellId; - - /** - * Type of reverse region-to-cell range bounds and - * iterators. - */ - typedef typename std::vector::const_iterator CellIter; - - /** - * Range of cells. Result from reverse (region-to-cell) - * mapping. - */ - class CellRange { - public: - /** - * Constructor. - * - * \param[in] b Beginning of range. - * \param[in] e One past end of range. - */ - CellRange(const CellIter b, - const CellIter e) - : b_(b), e_(e) - {} - - /** - * Read-only iterator on cell ranges. - */ - typedef CellIter const_iterator; - - /** - * Beginning of cell range. - */ - const_iterator begin() const { return b_; } - - /** - * One past end of cell range. - */ - const_iterator end() const { return e_; } - - private: - const_iterator b_; - const_iterator e_; - }; - - /** - * Number of declared regions in cell-to-region mapping. - */ - RegionId - numRegions() const { return RegionId(rev_.p.size()) - 1; } - - /** - * Compute region number of given active cell. - * - * \param[in] c Active cell - * \return Region to which @c c belongs. - */ - RegionId - region(const CellId c) const { return reg_[c]; } - - /** - * Extract active cells in particular region. - * - * \param[in] r Region number - * \returns Range of active cells in region @c r. - */ - CellRange - cells(const RegionId r) const { - const RegionId i = r - rev_.low; - return CellRange(rev_.c.begin() + rev_.p[i + 0], - rev_.c.begin() + rev_.p[i + 1]); - } - - private: - /** - * Copy of forward region mapping (cell-to-region). - */ - Region reg_; - - /** - * Reverse mapping (region-to-cell). - */ - struct { - typedef typename std::vector::size_type Pos; - std::vector p; /**< Region start pointers */ - std::vector c; /**< Region cells */ - RegionId low; /**< Smallest region number */ - - /** - * Compute reverse mapping. Standard linear insertion - * sort algorithm. - */ - void - init(const Region& reg) - { - typedef typename Region::const_iterator CI; - const std::pair - m = std::minmax_element(reg.begin(), reg.end()); - - low = *m.first; - - const typename Region::size_type - n = *m.second - low + 1; - - p.resize(n + 1); std::fill(p.begin(), p.end(), Pos(0)); - for (CellId i = 0, nc = reg.size(); i < nc; ++i) { - p[ reg[i] - low + 1 ] += 1; - } - - for (typename std::vector::size_type - i = 1, sz = p.size(); i < sz; ++i) { - p[0] += p[i]; - p[i] = p[0] - p[i]; - } - - assert (p[0] == - static_cast(reg.size())); - - c.resize(reg.size()); - for (CellId i = 0, nc = reg.size(); i < nc; ++i) { - c[ p[ reg[i] - low + 1 ] ++ ] = i; - } - - p[0] = 0; - } - } rev_; /**< Reverse mapping instance */ - }; /** * Equilibration record. diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index c4db1786a..4519f5f58 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -287,12 +287,12 @@ BOOST_AUTO_TEST_CASE (RegMapping) G->cartdims, cdim, &cells[0], &eqlnum[0]); } - Opm::equil::RegionMapping<> eqlmap(eqlnum); + Opm::RegionMapping<> eqlmap(eqlnum); PPress ppress(2, PVal(G->number_of_cells, 0)); for (int r = 0, e = eqlmap.numRegions(); r != e; ++r) { - const Opm::equil::RegionMapping<>::CellRange& + const Opm::RegionMapping<>::CellRange& rng = eqlmap.cells(r); const int rno = r; @@ -301,7 +301,7 @@ BOOST_AUTO_TEST_CASE (RegMapping) Opm::equil::phasePressures(*G, region[rno], rng, grav); PVal::size_type i = 0; - for (Opm::equil::RegionMapping<>::CellRange::const_iterator + for (Opm::RegionMapping<>::CellRange::const_iterator c = rng.begin(), ce = rng.end(); c != ce; ++c, ++i) { From 64b6a40191b4c64c65e8de490b5195fdcd84912a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 24 Feb 2014 15:55:14 +0100 Subject: [PATCH 29/53] Capitalize nested namespace names. equil -> Equil miscibility -> Miscibility --- opm/core/simulator/initStateEquil.hpp | 16 ++-- opm/core/simulator/initStateEquil_impl.hpp | 4 +- tests/test_equil.cpp | 88 +++++++++++----------- 3 files changed, 54 insertions(+), 54 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index fd1c6721a..1f6f11c5d 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -49,7 +49,7 @@ namespace Opm * This namespace is intentionally nested to avoid name clashes * with other parts of OPM. */ - namespace equil { + namespace Equil { template class DensityCalculator; @@ -116,7 +116,7 @@ namespace Opm * Types and routines relating to phase mixing in * equilibration calculations. */ - namespace miscibility { + namespace Miscibility { /** * Type that implements "no phase mixing" policy. */ @@ -258,7 +258,7 @@ namespace Opm return A_[np*opos + gpos] / A_[np*opos + opos]; } }; - } // namespace miscibility + } // namespace Miscibility @@ -336,8 +336,8 @@ namespace Opm * depth and (gas) pressure @c press. */ template + class RS = Miscibility::NoMixing, + class RV = Miscibility::NoMixing> class EquilReg { public: /** @@ -825,7 +825,7 @@ namespace Opm const UnstructuredGrid& G , const double grav) { - typedef miscibility::NoMixing NoMix; + typedef Miscibility::NoMixing NoMix; for (typename RMap::RegionId r = 0, nr = reg.numRegions(); @@ -863,7 +863,7 @@ namespace Opm const UnstructuredGrid& G , const double grav) { - typedef miscibility::NoMixing NoMix; + typedef Miscibility::NoMixing NoMix; for (typename RMap::RegionId r = 0, nr = reg.numRegions(); @@ -908,7 +908,7 @@ namespace Opm }; } // namespace DeckDependent - } // namespace equil + } // namespace Equil } // namespace Opm #include diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index a713e3b34..29ab6ac1d 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -485,7 +485,7 @@ namespace Opm } } // namespace Details - namespace equil { + namespace Equil { template std::vector< std::vector > @@ -568,7 +568,7 @@ namespace Opm return press; } - } // namespace equil + } // namespace Equil } // namespace Opm #endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 4519f5f58..88bec0729 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -64,27 +64,27 @@ BOOST_AUTO_TEST_CASE (PhasePressure) typedef Opm::BlackoilPropertiesBasic Props; Props props(param, G->dimensions, G->number_of_cells); - typedef Opm::equil::DensityCalculator RhoCalc; + typedef Opm::Equil::DensityCalculator RhoCalc; RhoCalc calc(props, 0); - Opm::equil::EquilRecord record = + Opm::Equil::EquilRecord record = { { 0 , 1e5 } , // Datum depth, pressure { 5 , 0 } , // Zwoc , Pcow_woc { 0 , 0 } // Zgoc , Pcgo_goc }; - Opm::equil::EquilReg + Opm::Equil::EquilReg region(record, calc, - Opm::equil::miscibility::NoMixing(), - Opm::equil::miscibility::NoMixing(), + Opm::Equil::Miscibility::NoMixing(), + Opm::Equil::Miscibility::NoMixing(), 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 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; @@ -118,10 +118,10 @@ BOOST_AUTO_TEST_CASE (CellSubset) typedef Opm::BlackoilPropertiesBasic Props; Props props(param, G->dimensions, G->number_of_cells); - typedef Opm::equil::DensityCalculator RhoCalc; + typedef Opm::Equil::DensityCalculator RhoCalc; RhoCalc calc(props, 0); - Opm::equil::EquilRecord record[] = + Opm::Equil::EquilRecord record[] = { { { 0 , 1e5 } , // Datum depth, pressure @@ -136,26 +136,26 @@ BOOST_AUTO_TEST_CASE (CellSubset) } }; - Opm::equil::EquilReg region[] = + Opm::Equil::EquilReg region[] = { - Opm::equil::EquilReg(record[0], calc, - Opm::equil::miscibility::NoMixing(), - Opm::equil::miscibility::NoMixing(), + Opm::Equil::EquilReg(record[0], calc, + Opm::Equil::Miscibility::NoMixing(), + Opm::Equil::Miscibility::NoMixing(), props.phaseUsage()) , - Opm::equil::EquilReg(record[0], calc, - Opm::equil::miscibility::NoMixing(), - Opm::equil::miscibility::NoMixing(), + Opm::Equil::EquilReg(record[0], calc, + Opm::Equil::Miscibility::NoMixing(), + Opm::Equil::Miscibility::NoMixing(), props.phaseUsage()) , - Opm::equil::EquilReg(record[1], calc, - Opm::equil::miscibility::NoMixing(), - Opm::equil::miscibility::NoMixing(), + Opm::Equil::EquilReg(record[1], calc, + Opm::Equil::Miscibility::NoMixing(), + Opm::Equil::Miscibility::NoMixing(), props.phaseUsage()) , - Opm::equil::EquilReg(record[1], calc, - Opm::equil::miscibility::NoMixing(), - Opm::equil::miscibility::NoMixing(), + Opm::Equil::EquilReg(record[1], calc, + Opm::Equil::Miscibility::NoMixing(), + Opm::Equil::Miscibility::NoMixing(), props.phaseUsage()) }; @@ -187,7 +187,7 @@ BOOST_AUTO_TEST_CASE (CellSubset) const int rno = int(r - cells.begin()); const double grav = 10; const PPress p = - Opm::equil::phasePressures(*G, region[rno], *r, grav); + Opm::Equil::phasePressures(*G, region[rno], *r, grav); PVal::size_type i = 0; for (std::vector::const_iterator @@ -233,10 +233,10 @@ BOOST_AUTO_TEST_CASE (RegMapping) typedef Opm::BlackoilPropertiesBasic Props; Props props(param, G->dimensions, G->number_of_cells); - typedef Opm::equil::DensityCalculator RhoCalc; + typedef Opm::Equil::DensityCalculator RhoCalc; RhoCalc calc(props, 0); - Opm::equil::EquilRecord record[] = + Opm::Equil::EquilRecord record[] = { { { 0 , 1e5 } , // Datum depth, pressure @@ -251,26 +251,26 @@ BOOST_AUTO_TEST_CASE (RegMapping) } }; - Opm::equil::EquilReg region[] = + Opm::Equil::EquilReg region[] = { - Opm::equil::EquilReg(record[0], calc, - Opm::equil::miscibility::NoMixing(), - Opm::equil::miscibility::NoMixing(), + Opm::Equil::EquilReg(record[0], calc, + Opm::Equil::Miscibility::NoMixing(), + Opm::Equil::Miscibility::NoMixing(), props.phaseUsage()) , - Opm::equil::EquilReg(record[0], calc, - Opm::equil::miscibility::NoMixing(), - Opm::equil::miscibility::NoMixing(), + Opm::Equil::EquilReg(record[0], calc, + Opm::Equil::Miscibility::NoMixing(), + Opm::Equil::Miscibility::NoMixing(), props.phaseUsage()) , - Opm::equil::EquilReg(record[1], calc, - Opm::equil::miscibility::NoMixing(), - Opm::equil::miscibility::NoMixing(), + Opm::Equil::EquilReg(record[1], calc, + Opm::Equil::Miscibility::NoMixing(), + Opm::Equil::Miscibility::NoMixing(), props.phaseUsage()) , - Opm::equil::EquilReg(record[1], calc, - Opm::equil::miscibility::NoMixing(), - Opm::equil::miscibility::NoMixing(), + Opm::Equil::EquilReg(record[1], calc, + Opm::Equil::Miscibility::NoMixing(), + Opm::Equil::Miscibility::NoMixing(), props.phaseUsage()) }; @@ -298,7 +298,7 @@ BOOST_AUTO_TEST_CASE (RegMapping) const int rno = r; const double grav = 10; const PPress p = - Opm::equil::phasePressures(*G, region[rno], rng, grav); + Opm::Equil::phasePressures(*G, region[rno], rng, grav); PVal::size_type i = 0; for (Opm::RegionMapping<>::CellRange::const_iterator @@ -328,7 +328,7 @@ BOOST_AUTO_TEST_CASE (DeckAllDead) grid(create_grid_cart3d(1, 1, 10), destroy_grid); Opm::EclipseGridParser deck("deadfluids.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, *grid, false); - Opm::equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, *grid, 10.0); + Opm::Equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, *grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells); @@ -364,7 +364,7 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion) const std::vector s = { 0.2, 0.2, 0.2, 0.466666666666, 0.733333333333, 1.0, 1.0, 1.0, 1.0 }; BOOST_REQUIRE(pc.size() == s.size()); for (size_t i = 0; i < pc.size(); ++i) { - const double s_computed = Opm::equil::satFromPc(props, phase, cell, pc[i], increasing); + const double s_computed = Opm::Equil::satFromPc(props, phase, cell, pc[i], increasing); BOOST_CHECK_CLOSE(s_computed, s[i], reltol); } } @@ -377,7 +377,7 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion) const std::vector s = { 0.8, 0.8, 0.8, 0.533333333333, 0.266666666666, 0.0, 0.0, 0.0, 0.0 }; BOOST_REQUIRE(pc.size() == s.size()); for (size_t i = 0; i < pc.size(); ++i) { - const double s_computed = Opm::equil::satFromPc(props, phase, cell, pc[i], increasing); + const double s_computed = Opm::Equil::satFromPc(props, phase, cell, pc[i], increasing); BOOST_CHECK_CLOSE(s_computed, s[i], reltol); } } @@ -390,7 +390,7 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion) const std::vector s = { 0.2, 0.333333333333, 0.6, 0.866666666666, 1.0 }; BOOST_REQUIRE(pc.size() == s.size()); for (size_t i = 0; i < pc.size(); ++i) { - const double s_computed = Opm::equil::satFromSumOfPcs(props, water, gas, cell, pc[i]); + const double s_computed = Opm::Equil::satFromSumOfPcs(props, water, gas, cell, pc[i]); BOOST_CHECK_CLOSE(s_computed, s[i], reltol); } } @@ -405,7 +405,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary) Opm::EclipseGridParser deck("capillary.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); - Opm::equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, grid, 10.0); + Opm::Equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -443,7 +443,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap) Opm::EclipseGridParser deck("capillary_overlap.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); - Opm::equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, grid, 10.0); + Opm::Equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); From 7e0bd6220523b8280b2d09896c4dab453eb02701 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 24 Feb 2014 16:09:04 +0100 Subject: [PATCH 30/53] Moved equilibration utilities to separate file. --- opm/core/simulator/EquilibrationHelpers.hpp | 631 ++++++++++++++++++++ opm/core/simulator/initStateEquil.hpp | 540 +---------------- 2 files changed, 632 insertions(+), 539 deletions(-) create mode 100644 opm/core/simulator/EquilibrationHelpers.hpp diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp new file mode 100644 index 000000000..442702d6c --- /dev/null +++ b/opm/core/simulator/EquilibrationHelpers.hpp @@ -0,0 +1,631 @@ +/* + 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 + + +/* +---- synopsis of EquilibrationHelpers.hpp ---- + +namespace Opm +{ + namespace Equil { + + template + class DensityCalculator; + + template <> + class DensityCalculator< BlackoilPropertiesInterface >; + + namespace Miscibility { + struct 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 { + /** + * Type that implements "no phase mixing" policy. + */ + struct NoMixing { + /** + * 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 + { + 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: + /** + * Constructor. + * + * \param[in] depth Depth nodes. + * \param[in] rs Dissolved gas-oil ratio at @c depth. + */ + RsVD(const std::vector& depth, + const std::vector& rs) + : depth_(depth) + , rs_(rs) + { + } + + /** + * 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 + { + return linearInterpolation(depth_, rs_, depth); + } + + private: + std::vector depth_; /**< Depth nodes */ + std::vector rs_; /**< Dissolved gas-oil ratio */ + }; + + + /** + * 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: + /** + * 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 + { + return std::max(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]; + } + }; + } // 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; + }; + + /** + * 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. + * + * \tparam RS Type that provides access to a calculator for + * (initial) dissolved gas-oil ratios as a function of depth + * and (oil) pressure. Must implement an operator() declared + * as + * + * double + * operator()(const double depth, + * const double press) + * + * that calculates the dissolved gas-oil ratio at depth @c + * depth and (oil) pressure @c press. + * + * \tparam RV Type that provides access to a calculator for + * (initial) vapourised oil-gas ratios as a function of depth + * and (gas) pressure. Must implement an operator() declared + * as + * + * double + * operator()(const double depth, + * const double press) + * + * that calculates the vapourised oil-gas ratio at depth @c + * depth and (gas) 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, + const RS& rs, + const RV& 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 RS CalcDissolution; + + /** + * Type of vapourised oil-gas ratio calculator. + */ + typedef RV 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 */ + RS rs_; /**< RS calculator */ + RV 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 index 1f6f11c5d..8e0490175 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -20,6 +20,7 @@ #ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED #define OPM_INITSTATEEQUIL_HEADER_INCLUDED +#include #include #include #include @@ -50,545 +51,6 @@ namespace Opm * 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 { - /** - * Type that implements "no phase mixing" policy. - */ - struct NoMixing { - /** - * 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 - { - 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: - /** - * Constructor. - * - * \param[in] depth Depth nodes. - * \param[in] rs Dissolved gas-oil ratio at @c depth. - */ - RsVD(const std::vector& depth, - const std::vector& rs) - : depth_(depth) - , rs_(rs) - { - } - - /** - * 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 - { - return linearInterpolation(depth_, rs_, depth); - } - - private: - std::vector depth_; /**< Depth nodes */ - std::vector rs_; /**< Dissolved gas-oil ratio */ - }; - - - /** - * 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: - /** - * 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 - { - return std::max(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]; - } - }; - } // 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; - }; - - /** - * 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. - * - * \tparam RS Type that provides access to a calculator for - * (initial) dissolved gas-oil ratios as a function of depth - * and (oil) pressure. Must implement an operator() declared - * as - * - * double - * operator()(const double depth, - * const double press) - * - * that calculates the dissolved gas-oil ratio at depth @c - * depth and (oil) pressure @c press. - * - * \tparam RV Type that provides access to a calculator for - * (initial) vapourised oil-gas ratios as a function of depth - * and (gas) pressure. Must implement an operator() declared - * as - * - * double - * operator()(const double depth, - * const double press) - * - * that calculates the vapourised oil-gas ratio at depth @c - * depth and (gas) 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, - const RS& rs, - const RV& 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 RS CalcDissolution; - - /** - * Type of vapourised oil-gas ratio calculator. - */ - typedef RV 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 */ - RS rs_; /**< RS calculator */ - RV 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; - } - } - - /** * Compute initial phase pressures by means of equilibration. From 74573947bb0953bc0843c485afa849203ddc256b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 24 Feb 2014 16:11:50 +0100 Subject: [PATCH 31/53] Prune includes. --- opm/core/simulator/initStateEquil.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 8e0490175..9be77685d 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -24,9 +24,7 @@ #include #include #include -#include #include -#include #include #include From 8266c52d996697252fa752943dcb04d5c25782c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Feb 2014 14:16:51 +0100 Subject: [PATCH 32/53] Made NoMixing a class. For uniformity with its sibling classes. --- opm/core/simulator/EquilibrationHelpers.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index 442702d6c..316b435be 100644 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ b/opm/core/simulator/EquilibrationHelpers.hpp @@ -41,7 +41,7 @@ namespace Opm class DensityCalculator< BlackoilPropertiesInterface >; namespace Miscibility { - struct NoMixing; + class NoMixing; class RsVD; class RsSatAtContact; } @@ -156,7 +156,8 @@ namespace Opm /** * Type that implements "no phase mixing" policy. */ - struct NoMixing { + class NoMixing { + public: /** * Function call. * From 38c89f363d21a4c1cd3a1d7175b01da812a09b91 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Feb 2014 14:47:24 +0100 Subject: [PATCH 33/53] Made phase mixing functors a class hierarchy. In summary: - added RsFunction (base class), - made NoMixing, RsVD, RsSatAtContact inherit RsFunction, - RS and RV are no longer template arguments for EquilReg class, - EquilReg constructor now takes two shared_ptr, - use of constructor updated, mostly using make_shared. --- opm/core/simulator/EquilibrationHelpers.hpp | 84 +++++++++++---------- opm/core/simulator/initStateEquil.hpp | 6 +- tests/test_equil.cpp | 36 ++++----- 3 files changed, 65 insertions(+), 61 deletions(-) diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index 316b435be..0588f3464 100644 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ b/opm/core/simulator/EquilibrationHelpers.hpp @@ -26,6 +26,8 @@ #include #include +#include + /* ---- synopsis of EquilibrationHelpers.hpp ---- @@ -41,6 +43,7 @@ namespace Opm class DensityCalculator< BlackoilPropertiesInterface >; namespace Miscibility { + class RsFunction; class NoMixing; class RsVD; class RsSatAtContact; @@ -48,9 +51,7 @@ namespace Opm struct EquilRecord; - template + template class EquilReg; @@ -148,15 +149,40 @@ namespace Opm 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 = 0; + }; + + /** * Type that implements "no phase mixing" policy. */ - class NoMixing { + class NoMixing : public RsFunction { public: /** * Function call. @@ -179,12 +205,13 @@ namespace Opm } }; + /** * Type that implements "dissolved gas-oil ratio" * tabulated as a function of depth policy. Data * typically taken from keyword 'RSVD'. */ - class RsVD { + class RsVD : public RsFunction { public: /** * Constructor. @@ -238,7 +265,7 @@ namespace Opm * This should yield Rs-values that are constant below the * contact, and decreasing above the contact. */ - class RsSatAtContact { + class RsSatAtContact : public RsFunction { public: /** * Constructor. @@ -295,6 +322,7 @@ namespace Opm return A_[np*opos + gpos] / A_[np*opos + opos]; } }; + } // namespace Miscibility @@ -347,34 +375,8 @@ namespace Opm * * that calculates the phase densities of all phases in @c * svol at fluid pressure @c press. - * - * \tparam RS Type that provides access to a calculator for - * (initial) dissolved gas-oil ratios as a function of depth - * and (oil) pressure. Must implement an operator() declared - * as - * - * double - * operator()(const double depth, - * const double press) - * - * that calculates the dissolved gas-oil ratio at depth @c - * depth and (oil) pressure @c press. - * - * \tparam RV Type that provides access to a calculator for - * (initial) vapourised oil-gas ratios as a function of depth - * and (gas) pressure. Must implement an operator() declared - * as - * - * double - * operator()(const double depth, - * const double press) - * - * that calculates the vapourised oil-gas ratio at depth @c - * depth and (gas) pressure @c press. */ - template + template class EquilReg { public: /** @@ -388,8 +390,8 @@ namespace Opm */ EquilReg(const EquilRecord& rec, const DensCalc& density, - const RS& rs, - const RV& rv, + std::shared_ptr rs, + std::shared_ptr rv, const PhaseUsage& pu) : rec_ (rec) , density_(density) @@ -407,12 +409,12 @@ namespace Opm /** * Type of dissolved gas-oil ratio calculator. */ - typedef RS CalcDissolution; + typedef Miscibility::RsFunction CalcDissolution; /** * Type of vapourised oil-gas ratio calculator. */ - typedef RV CalcEvaporation; + typedef Miscibility::RsFunction CalcEvaporation; /** * Datum depth in current region @@ -459,14 +461,14 @@ namespace Opm * region. */ const CalcDissolution& - dissolutionCalculator() const { return this->rs_; } + dissolutionCalculator() const { return *this->rs_; } /** * Retrieve vapourised oil-gas ratio calculator of current * region. */ const CalcEvaporation& - evaporationCalculator() const { return this->rv_; } + evaporationCalculator() const { return *this->rv_; } /** * Retrieve active fluid phase summary. @@ -477,8 +479,8 @@ namespace Opm private: EquilRecord rec_; /**< Equilibration data */ DensCalc density_; /**< Density calculator */ - RS rs_; /**< RS calculator */ - RV rv_; /**< RV calculator */ + std::shared_ptr rs_; /**< RS calculator */ + std::shared_ptr rv_; /**< RV calculator */ PhaseUsage pu_; /**< Active phase summary */ }; diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 9be77685d..aa9f78625 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -296,7 +296,8 @@ namespace Opm const int repcell = *cells.begin(); const RhoCalc calc(props, repcell); - const EqReg eqreg(rec[r], calc, NoMix(), NoMix(), + const EqReg eqreg(rec[r], calc, + std::make_shared(), std::make_shared(), props.phaseUsage()); const PPress& res = phasePressures(G, eqreg, cells, grav); @@ -334,7 +335,8 @@ namespace Opm const int repcell = *cells.begin(); const RhoCalc calc(props, repcell); - const EqReg eqreg(rec[r], calc, NoMix(), NoMix(), + const EqReg eqreg(rec[r], calc, + std::make_shared(), std::make_shared(), props.phaseUsage()); const PPress press = phasePressures(G, eqreg, cells, grav); diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 88bec0729..5a5e9019b 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -76,8 +76,8 @@ BOOST_AUTO_TEST_CASE (PhasePressure) Opm::Equil::EquilReg region(record, calc, - Opm::Equil::Miscibility::NoMixing(), - Opm::Equil::Miscibility::NoMixing(), + std::make_shared(), + std::make_shared(), props.phaseUsage()); std::vector cells(G->number_of_cells); @@ -139,23 +139,23 @@ BOOST_AUTO_TEST_CASE (CellSubset) Opm::Equil::EquilReg region[] = { Opm::Equil::EquilReg(record[0], calc, - Opm::Equil::Miscibility::NoMixing(), - Opm::Equil::Miscibility::NoMixing(), + std::make_shared(), + std::make_shared(), props.phaseUsage()) , Opm::Equil::EquilReg(record[0], calc, - Opm::Equil::Miscibility::NoMixing(), - Opm::Equil::Miscibility::NoMixing(), + std::make_shared(), + std::make_shared(), props.phaseUsage()) , Opm::Equil::EquilReg(record[1], calc, - Opm::Equil::Miscibility::NoMixing(), - Opm::Equil::Miscibility::NoMixing(), + std::make_shared(), + std::make_shared(), props.phaseUsage()) , Opm::Equil::EquilReg(record[1], calc, - Opm::Equil::Miscibility::NoMixing(), - Opm::Equil::Miscibility::NoMixing(), + std::make_shared(), + std::make_shared(), props.phaseUsage()) }; @@ -254,23 +254,23 @@ BOOST_AUTO_TEST_CASE (RegMapping) Opm::Equil::EquilReg region[] = { Opm::Equil::EquilReg(record[0], calc, - Opm::Equil::Miscibility::NoMixing(), - Opm::Equil::Miscibility::NoMixing(), + std::make_shared(), + std::make_shared(), props.phaseUsage()) , Opm::Equil::EquilReg(record[0], calc, - Opm::Equil::Miscibility::NoMixing(), - Opm::Equil::Miscibility::NoMixing(), + std::make_shared(), + std::make_shared(), props.phaseUsage()) , Opm::Equil::EquilReg(record[1], calc, - Opm::Equil::Miscibility::NoMixing(), - Opm::Equil::Miscibility::NoMixing(), + std::make_shared(), + std::make_shared(), props.phaseUsage()) , Opm::Equil::EquilReg(record[1], calc, - Opm::Equil::Miscibility::NoMixing(), - Opm::Equil::Miscibility::NoMixing(), + std::make_shared(), + std::make_shared(), props.phaseUsage()) }; From d5ffcc051c2706edbe1f88b14010fef73496868c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Feb 2014 14:49:06 +0100 Subject: [PATCH 34/53] Removed redundant calcPressII() method. Pressure is also calculated in the calcPressSat() method. --- opm/core/simulator/initStateEquil.hpp | 53 +++------------------------ 1 file changed, 6 insertions(+), 47 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index aa9f78625..9ba625244 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -259,9 +259,7 @@ namespace Opm { const std::vector rec = getEquil(deck); const RegionMapping<> eqlmap(equilnum(deck, G)); - - calcPressII(eqlmap, rec, props, G, grav); - calcSat(eqlmap, rec, props, G, grav); + calcPressSat(eqlmap, rec, props, G, grav); } typedef std::vector PVal; @@ -279,50 +277,11 @@ namespace Opm template void - calcPressII(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, - std::make_shared(), std::make_shared(), - props.phaseUsage()); - - const PPress& res = phasePressures(G, eqreg, cells, grav); - - for (int p = 0, np = props.numPhases(); p < np; ++p) { - PVal& d = pp_[p]; - PVal::const_iterator s = res[p].begin(); - for (typename RMap::CellRange::const_iterator - c = cells.begin(), - e = cells.end(); - c != e; ++c, ++s) - { - d[*c] = *s; - } - } - } - } - - template - void - calcSat(const RMap& reg , - const std::vector< EquilRecord >& rec , - const Opm::BlackoilPropertiesInterface& props, - const UnstructuredGrid& G , - const double grav) + calcPressSat(const RMap& reg , + const std::vector< EquilRecord >& rec , + const Opm::BlackoilPropertiesInterface& props, + const UnstructuredGrid& G , + const double grav) { typedef Miscibility::NoMixing NoMix; From 27743b4a1387835a3c78484474235e355a4da0a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 08:31:03 +0100 Subject: [PATCH 35/53] Enable live oil in initialisation. --- opm/core/simulator/initStateEquil.hpp | 29 ++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 9ba625244..d2ebae79f 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -257,8 +257,33 @@ namespace Opm sat_(props.numPhases(), std::vector(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")) { + if (deck.hasField("RSVD")) { + // Rs has been specified as a function of depth. + OPM_THROW(std::runtime_error, "Cannot initialise: RSVD field not read by EclipseGridParser class."); + } else { + // Default initialisation: constant Rs below contact, saturated above. + for (size_t i = 0; i < rec.size(); ++i) { + const int cell = *(eqlmap.cells(i + 1).begin()); + const double p_contact = rec[i].goc.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()); + } + } + + // Compute phase pressures and saturations. calcPressSat(eqlmap, rec, props, G, grav); } @@ -272,6 +297,8 @@ namespace Opm typedef DensityCalculator RhoCalc; typedef EquilReg EqReg; + std::vector< std::shared_ptr > rs_func_; + PPress pp_; PPress sat_; @@ -295,7 +322,7 @@ namespace Opm const RhoCalc calc(props, repcell); const EqReg eqreg(rec[r], calc, - std::make_shared(), std::make_shared(), + rs_func_[r], std::make_shared(), props.phaseUsage()); const PPress press = phasePressures(G, eqreg, cells, grav); From 91ae9df907870f2e8f66b81d2b76eb53b9340f8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 08:31:51 +0100 Subject: [PATCH 36/53] Add test for live oil initialisation. The test is not finished or verified yet. --- tests/equil_liveoil.DATA | 46 +++++++++++++++++++++++++++++++++++++++ tests/test_equil.cpp | 47 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 tests/equil_liveoil.DATA diff --git a/tests/equil_liveoil.DATA b/tests/equil_liveoil.DATA new file mode 100644 index 000000000..5b2766915 --- /dev/null +++ b/tests/equil_liveoil.DATA @@ -0,0 +1,46 @@ +WATER +OIL +GAS +DISGAS + + + +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 +50 150 50 0.25 45 0.35 +/ diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 5a5e9019b..3e1866772 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -480,4 +480,51 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap) } + +BOOST_AUTO_TEST_CASE (DeckWithLiveOil) +{ + Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0); + const UnstructuredGrid& grid = *(gm.c_grid()); + Opm::EclipseGridParser deck("equil_liveoil.DATA"); + Opm::BlackoilPropertiesFromDeck props(deck, grid, false); + + Opm::Equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, grid, 10.0); + const auto& pressures = comp.press(); + BOOST_REQUIRE(pressures.size() == 3); + BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-6; + BOOST_CHECK_CLOSE(pressures[0][first] , 1.45e7 , reltol); + BOOST_CHECK_CLOSE(pressures[0][last ] , 1.545e7 , reltol); + BOOST_CHECK_CLOSE(pressures[1][last] , 1.5489764605846416e7 , reltol); + + const auto& sats = comp.saturation(); + // std::cout << "Saturations:\n"; + // for (const auto& sat : sats) { + // for (const double s : sat) { + // std::cout << s << ' '; + // } + // std::cout << std::endl; + // } + const std::vector s[3]{ + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.223141818182, 0.532269090909, 0.78471, 0.91526, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.207743333333, 0.08474, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.776858181818, 0.467730909091, 0.0075466666666, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + for (int phase = 0; phase < 3; ++phase) { + BOOST_REQUIRE(sats[phase].size() == s[phase].size()); + for (size_t i = 0; i < s[phase].size(); ++i) { + std::cout << sats[phase][i] << '\n'; + //BOOST_CHECK_CLOSE(sats[phase][i], s[phase][i], reltol); + } + std::cout << std::endl; + } +} + + BOOST_AUTO_TEST_SUITE_END() From a7c45d4e9ff9cf894045e7735d373b3e7aa1b4ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 09:08:39 +0100 Subject: [PATCH 37/53] Rename PhasePressureSaturationComputer -> InitialStateComputer. Also add (unused so far) rs_ field to class. --- opm/core/simulator/initStateEquil.hpp | 51 +++++++++++++++++---------- tests/test_equil.cpp | 8 ++--- 2 files changed, 36 insertions(+), 23 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index d2ebae79f..ee8bf7ac6 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -243,19 +243,20 @@ namespace Opm } template - class PhasePressureSaturationComputer; + class InitialStateComputer; template <> - class PhasePressureSaturationComputer { + class InitialStateComputer { public: - PhasePressureSaturationComputer(const BlackoilPropertiesInterface& props, - const EclipseGridParser& deck , - const UnstructuredGrid& G , - const double grav = unit::gravity) + 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)) + std::vector(G.number_of_cells)), + rs_(G.number_of_cells) { // Get the equilibration records. const std::vector rec = getEquil(deck); @@ -287,11 +288,11 @@ namespace Opm calcPressSat(eqlmap, rec, props, G, grav); } - typedef std::vector PVal; - typedef std::vector PPress; + typedef std::vector Vec; + typedef std::vector PVec; // One per phase. - const PPress& press() const { return pp_; } - const PPress& saturation() const { return sat_; } + const PVec& press() const { return pp_; } + const PVec& saturation() const { return sat_; } private: typedef DensityCalculator RhoCalc; @@ -299,8 +300,9 @@ namespace Opm std::vector< std::shared_ptr > rs_func_; - PPress pp_; - PPress sat_; + PVec pp_; + PVec sat_; + Vec rs_; template void @@ -325,12 +327,13 @@ namespace Opm rs_func_[r], std::make_shared(), props.phaseUsage()); - const PPress press = phasePressures(G, eqreg, cells, grav); - const PPress sat = phaseSaturations(eqreg, cells, props, press); + const PVec press = phasePressures(G, eqreg, cells, grav); + const PVec sat = phaseSaturations(eqreg, cells, props, press); + const Vec rs(cells.size());// = gasOilRatio(); for (int p = 0, np = props.numPhases(); p < np; ++p) { - PVal& d = pp_[p]; - PVal::const_iterator s = press[p].begin(); + Vec& d = pp_[p]; + Vec::const_iterator s = press[p].begin(); for (typename RMap::CellRange::const_iterator c = cells.begin(), e = cells.end(); @@ -340,8 +343,8 @@ namespace Opm } } for (int p = 0, np = props.numPhases(); p < np; ++p) { - PVal& d = sat_[p]; - PVal::const_iterator s = sat[p].begin(); + Vec& d = sat_[p]; + Vec::const_iterator s = sat[p].begin(); for (typename RMap::CellRange::const_iterator c = cells.begin(), e = cells.end(); @@ -350,6 +353,16 @@ namespace Opm d[*c] = *s; } } + Vec::const_iterator s = rs.begin(); + Vec& d = rs_; + for (typename RMap::CellRange::const_iterator + c = cells.begin(), + e = cells.end(); + c != e; ++c, ++s) + { + d[*c] = *s; + } + } } diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 3e1866772..38bdff80d 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -328,7 +328,7 @@ BOOST_AUTO_TEST_CASE (DeckAllDead) grid(create_grid_cart3d(1, 1, 10), destroy_grid); Opm::EclipseGridParser deck("deadfluids.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, *grid, false); - Opm::Equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, *grid, 10.0); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, *grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells); @@ -405,7 +405,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary) Opm::EclipseGridParser deck("capillary.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); - Opm::Equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, grid, 10.0); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -443,7 +443,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap) Opm::EclipseGridParser deck("capillary_overlap.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); - Opm::Equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, grid, 10.0); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -488,7 +488,7 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil) Opm::EclipseGridParser deck("equil_liveoil.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); - Opm::Equil::DeckDependent::PhasePressureSaturationComputer comp(props, deck, grid, 10.0); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); From a1192c09dc7327f919cfbc0fea95ad86edb1d883 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 09:31:48 +0100 Subject: [PATCH 38/53] Refactor copying of region to global data. --- opm/core/simulator/initStateEquil.hpp | 48 ++++++++++----------------- 1 file changed, 17 insertions(+), 31 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index ee8bf7ac6..6ddb80ef8 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -331,41 +331,27 @@ namespace Opm const PVec sat = phaseSaturations(eqreg, cells, props, press); const Vec rs(cells.size());// = gasOilRatio(); - for (int p = 0, np = props.numPhases(); p < np; ++p) { - Vec& d = pp_[p]; - Vec::const_iterator s = press[p].begin(); - for (typename RMap::CellRange::const_iterator - c = cells.begin(), - e = cells.end(); - c != e; ++c, ++s) - { - d[*c] = *s; - } + const int np = props.numPhases(); + for (int p = 0; p < np; ++p) { + copyFromRegion(press[p], cells, pp_[p]); + copyFromRegion(sat[p], cells, sat_[p]); } - for (int p = 0, np = props.numPhases(); p < np; ++p) { - Vec& d = sat_[p]; - Vec::const_iterator s = sat[p].begin(); - for (typename RMap::CellRange::const_iterator - c = cells.begin(), - e = cells.end(); - c != e; ++c, ++s) - { - d[*c] = *s; - } - } - Vec::const_iterator s = rs.begin(); - Vec& d = rs_; - for (typename RMap::CellRange::const_iterator - c = cells.begin(), - e = cells.end(); - c != e; ++c, ++s) - { - d[*c] = *s; - } - + copyFromRegion(rs, cells, rs_); } } + 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 From 0fbc3a1ccb57a67d07519281048eecfd3c640218 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 09:37:48 +0100 Subject: [PATCH 39/53] Moved implementation of phaseSaturations() to _impl file. --- opm/core/simulator/initStateEquil.hpp | 58 +----------------- opm/core/simulator/initStateEquil_impl.hpp | 70 ++++++++++++++++++++++ 2 files changed, 71 insertions(+), 57 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 6ddb80ef8..296e67747 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -126,63 +126,7 @@ namespace Opm phaseSaturations(const Region& reg, const CellRange& cells, const BlackoilPropertiesInterface& props, - const 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; - } - 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; - } - phase_saturations[oilpos][local_index] = 1.0 - sw - sg; - } - return phase_saturations; - } - + const std::vector< std::vector >& phase_pressures); namespace DeckDependent { diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index 29ab6ac1d..03c10b4a4 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -485,7 +485,10 @@ namespace Opm } } // namespace Details + namespace Equil { + + template std::vector< std::vector > @@ -568,6 +571,73 @@ namespace Opm return press; } + + + + + template + std::vector< std::vector > + phaseSaturations(const Region& reg, + const CellRange& cells, + const BlackoilPropertiesInterface& props, + const 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; + } + 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; + } + phase_saturations[oilpos][local_index] = 1.0 - sw - sg; + } + return phase_saturations; + } + + } // namespace Equil } // namespace Opm From 4aaeff0078c6fbc79b8545a623713e78abf0a070 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 10:39:18 +0100 Subject: [PATCH 40/53] Added Rv field to InitialStateComputer. It is currently not computed, as for Rs. --- opm/core/simulator/initStateEquil.hpp | 37 ++++++++++++++++++++------- 1 file changed, 28 insertions(+), 9 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 296e67747..98842fb4b 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -129,7 +129,11 @@ namespace Opm const std::vector< std::vector >& phase_pressures); + + + namespace DeckDependent { + inline std::vector getEquil(const EclipseGridParser& deck) @@ -168,6 +172,9 @@ namespace Opm } } + + + inline std::vector equilnum(const EclipseGridParser& deck, @@ -186,6 +193,9 @@ namespace Opm return eqlnum; } + + + template class InitialStateComputer; @@ -200,7 +210,8 @@ namespace Opm std::vector(G.number_of_cells)), sat_(props.numPhases(), std::vector(G.number_of_cells)), - rs_(G.number_of_cells) + rs_(G.number_of_cells), + rv_(G.number_of_cells) { // Get the equilibration records. const std::vector rec = getEquil(deck); @@ -228,8 +239,11 @@ namespace Opm } } - // Compute phase pressures and saturations. - calcPressSat(eqlmap, rec, props, G, grav); + // 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; @@ -237,6 +251,8 @@ namespace Opm 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; @@ -247,14 +263,15 @@ namespace Opm PVec pp_; PVec sat_; Vec rs_; + Vec rv_; template void - calcPressSat(const RMap& reg , - const std::vector< EquilRecord >& rec , - const Opm::BlackoilPropertiesInterface& props, - const UnstructuredGrid& G , - const double grav) + calcPressSatRsRv(const RMap& reg , + const std::vector< EquilRecord >& rec , + const Opm::BlackoilPropertiesInterface& props, + const UnstructuredGrid& G , + const double grav) { typedef Miscibility::NoMixing NoMix; @@ -273,7 +290,8 @@ namespace Opm const PVec press = phasePressures(G, eqreg, cells, grav); const PVec sat = phaseSaturations(eqreg, cells, props, press); - const Vec rs(cells.size());// = gasOilRatio(); + const Vec rs(cells.size()); + const Vec rv(cells.size()); const int np = props.numPhases(); for (int p = 0; p < np; ++p) { @@ -281,6 +299,7 @@ namespace Opm copyFromRegion(sat[p], cells, sat_[p]); } copyFromRegion(rs, cells, rs_); + copyFromRegion(rv, cells, rv_); } } From f84162bc39dd308de203e5c82f7d05cad134297e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 10:40:14 +0100 Subject: [PATCH 41/53] Add initStateEquil() function. It is not quite complete yet for the following reasons: - it does not compute state.surfacevol(), - the InitialStateComputer class does not compute Rs or Rv, - it has not been verified. --- opm/core/simulator/initStateEquil.hpp | 24 +++++++++ opm/core/simulator/initStateEquil_impl.hpp | 58 ++++++++++++++++++++++ 2 files changed, 82 insertions(+) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 98842fb4b..43b34d2c7 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -21,6 +21,7 @@ #define OPM_INITSTATEEQUIL_HEADER_INCLUDED #include +#include #include #include #include @@ -41,6 +42,29 @@ 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); + + /** * Types and routines that collectively implement a basic * ECLIPSE-style equilibration-based initialisation scheme. diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index 03c10b4a4..2b2f99dba 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -639,6 +639,64 @@ namespace Opm } // 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. + } + + + } // namespace Opm #endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED From af9c1fc4cc3de6943d5d99b82dcce79ee0ca7918 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 13:14:48 +0100 Subject: [PATCH 42/53] Add computeRs() function and use from InitialStateComputer. --- opm/core/simulator/initStateEquil.hpp | 36 +++++++++++++++++++--- opm/core/simulator/initStateEquil_impl.hpp | 35 +++++++++++++++++++++ 2 files changed, 66 insertions(+), 5 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 43b34d2c7..cb4051e05 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -154,6 +154,28 @@ namespace Opm + /** + * 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); namespace DeckDependent { @@ -314,16 +336,20 @@ namespace Opm const PVec press = phasePressures(G, eqreg, cells, grav); const PVec sat = phaseSaturations(eqreg, cells, props, press); - const Vec rs(cells.size()); - const Vec rv(cells.size()); - const int np = props.numPhases(); for (int p = 0; p < np; ++p) { copyFromRegion(press[p], cells, pp_[p]); copyFromRegion(sat[p], cells, sat_[p]); } - copyFromRegion(rs, cells, rs_); - copyFromRegion(rv, cells, rv_); + if (props.phaseUsage().phase_used[BlackoilPhases::Liquid] + && props.phaseUsage().phase_used[BlackoilPhases::Vapour]) { + const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid]; + const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r])); + const Vec rv(cells.size(), 0.0); + std::cout << "rs.size() = " << rs.size() << std::endl; + copyFromRegion(rs, cells, rs_); + copyFromRegion(rv, cells, rv_); + } } } diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index 2b2f99dba..c075f4ecb 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -638,6 +638,41 @@ namespace Opm } + /** + * 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) + { + 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]); + } + return rs; + } + + + } // namespace Equil From 1b3cac2433471f134a8f9b7cb9026cc81c7f47c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 13:27:07 +0100 Subject: [PATCH 43/53] Remove debugging output. --- opm/core/simulator/initStateEquil.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index cb4051e05..1b81b07bd 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -346,7 +346,6 @@ namespace Opm const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid]; const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r])); const Vec rv(cells.size(), 0.0); - std::cout << "rs.size() = " << rs.size() << std::endl; copyFromRegion(rs, cells, rs_); copyFromRegion(rv, cells, rv_); } From 200b5372dd3ca584f36d88aa220ff2094ca5358f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 13:55:15 +0100 Subject: [PATCH 44/53] Added simple program example for initialisation. --- examples/compute_initial_state.cpp | 104 +++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 examples/compute_initial_state.cpp diff --git a/examples/compute_initial_state.cpp b/examples/compute_initial_state.cpp new file mode 100644 index 000000000..0b11ecf1c --- /dev/null +++ b/examples/compute_initial_state.cpp @@ -0,0 +1,104 @@ +/* + 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 + +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"); + 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; +} From 1c9675605cc540e1a1f74ceeb50de673bfc12b76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 14:48:14 +0100 Subject: [PATCH 45/53] Fix bug in RS initialisation. Also throw if default init is specified and datum != goc depth. --- opm/core/simulator/initStateEquil.hpp | 8 +++++++- tests/equil_liveoil.DATA | 2 +- tests/test_equil.cpp | 6 +++--- 3 files changed, 11 insertions(+), 5 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 1b81b07bd..3767563f3 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -275,7 +275,13 @@ namespace Opm // Default initialisation: constant Rs below contact, saturated above. for (size_t i = 0; i < rec.size(); ++i) { const int cell = *(eqlmap.cells(i + 1).begin()); - const double p_contact = rec[i].goc.press; + 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)); } } diff --git a/tests/equil_liveoil.DATA b/tests/equil_liveoil.DATA index 5b2766915..955a09b7f 100644 --- a/tests/equil_liveoil.DATA +++ b/tests/equil_liveoil.DATA @@ -42,5 +42,5 @@ DENSITY / EQUIL -50 150 50 0.25 45 0.35 +45 150 50 0.25 45 0.35 / diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 38bdff80d..6a175acc2 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -499,9 +499,9 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil) // solver, and it is unclear if we should check it against // the true answer or something else. const double reltol = 1.0e-6; - BOOST_CHECK_CLOSE(pressures[0][first] , 1.45e7 , reltol); - BOOST_CHECK_CLOSE(pressures[0][last ] , 1.545e7 , reltol); - BOOST_CHECK_CLOSE(pressures[1][last] , 1.5489764605846416e7 , reltol); + BOOST_CHECK_CLOSE(pressures[0][first] , 1.4551328443106037e7 , reltol); + BOOST_CHECK_CLOSE(pressures[0][last ] , 1.5501328443106037e7 , reltol); + BOOST_CHECK_CLOSE(pressures[1][last] , 1.5541598197355453e7 , reltol); const auto& sats = comp.saturation(); // std::cout << "Saturations:\n"; From d13bf03e1f32088dafb1d34186ce079848ea992c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 14:57:38 +0100 Subject: [PATCH 46/53] Bugfix in RsSatAtContact: use min(), not max(). Also modified test to match output. --- opm/core/simulator/EquilibrationHelpers.hpp | 2 +- tests/test_equil.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index 0588f3464..0d4d08196 100644 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ b/opm/core/simulator/EquilibrationHelpers.hpp @@ -300,7 +300,7 @@ namespace Opm operator()(const double /* depth */, const double press) const { - return std::max(satRs(press), rs_sat_contact_); + return std::min(satRs(press), rs_sat_contact_); } private: diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 6a175acc2..b4df1fae8 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -499,9 +499,9 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil) // solver, and it is unclear if we should check it against // the true answer or something else. const double reltol = 1.0e-6; - BOOST_CHECK_CLOSE(pressures[0][first] , 1.4551328443106037e7 , reltol); - BOOST_CHECK_CLOSE(pressures[0][last ] , 1.5501328443106037e7 , reltol); - BOOST_CHECK_CLOSE(pressures[1][last] , 1.5541598197355453e7 , reltol); + BOOST_CHECK_CLOSE(pressures[0][first], 1.4551302072306179e7, reltol); + BOOST_CHECK_CLOSE(pressures[0][last], 1.5501302072306179e7, reltol); + BOOST_CHECK_CLOSE(pressures[1][last], 1.5538684664272346e7, reltol); const auto& sats = comp.saturation(); // std::cout << "Saturations:\n"; From c2b647d2459f3685838cd6acb0347d551081af72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Feb 2014 15:55:08 +0100 Subject: [PATCH 47/53] Add test data file for compute_initial_state.cpp. --- tests/equil_liveoil_grid.DATA | 65 +++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 tests/equil_liveoil_grid.DATA 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 +/ From e1d2fa2088b3953d68103f2254b87b6ce58fb57c Mon Sep 17 00:00:00 2001 From: osae Date: Wed, 26 Mar 2014 14:08:39 +0100 Subject: [PATCH 48/53] Some adjustments to equil initialisation. - Saturations, phase pressures, and standard initialsation of RS and RV now agree to baseline. - Tables of RS and RV versus vertical depth (kw RSVD RVVD) have been hardcoded for testing (need new parser) and the calculations agree to baseline in the gas and oil zones. In the water zone there is some differences: Our code computes saturated RS and RV using the final phase pressures (these are modified to be consistent with saturations and capillary pressures) while the baseline uses unmodified phase pressures. --- opm/core/simulator/EquilibrationHelpers.hpp | 213 +++++++++++++++++++- opm/core/simulator/initStateEquil.hpp | 91 +++++++-- opm/core/simulator/initStateEquil_impl.hpp | 34 +++- 3 files changed, 306 insertions(+), 32 deletions(-) diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index 0d4d08196..d89dd8c3c 100644 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ b/opm/core/simulator/EquilibrationHelpers.hpp @@ -175,7 +175,8 @@ namespace Opm * depth and pressure @c press. */ virtual double operator()(const double depth, - const double press) const = 0; + const double press, + const double sat = 0.0) const = 0; }; @@ -199,7 +200,8 @@ namespace Opm */ double operator()(const double /* depth */, - const double /* press */) const + const double /* press */, + const double sat = 0.0) const { return 0.0; } @@ -216,14 +218,24 @@ namespace Opm /** * 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 std::vector& depth, + RsVD(const BlackoilPropertiesInterface& props, + const int cell, + const std::vector& depth, const std::vector& rs) - : depth_(depth) + : 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; } /** @@ -240,14 +252,111 @@ namespace Opm */ double operator()(const double depth, - const double /* press */) const + const double press, + const double sat_gas = 0.0) const { - return linearInterpolation(depth_, rs_, depth); + 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]; + } }; @@ -298,9 +407,14 @@ namespace Opm */ double operator()(const double /* depth */, - const double press) const - { - return std::min(satRs(press), rs_sat_contact_); + 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: @@ -323,6 +437,84 @@ namespace Opm } }; + + /** + * 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 @@ -355,6 +547,9 @@ namespace Opm double depth; double press; } main, woc, goc; + int live_oil_table_index; + int wet_gas_table_index; + int N; }; /** diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 3767563f3..566f2b00f 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -175,8 +175,8 @@ namespace Opm std::vector computeRs(const UnstructuredGrid& grid, const CellRangeType& cells, const std::vector oil_pressure, - const Miscibility::RsFunction& rs_func); - + const Miscibility::RsFunction& rs_func, + const std::vector gas_saturation); namespace DeckDependent { @@ -205,8 +205,17 @@ namespace Opm , { 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); } @@ -267,14 +276,25 @@ namespace Opm // Create Rs functions. rs_func_.reserve(rec.size()); - if (deck.hasField("DISGAS")) { - if (deck.hasField("RSVD")) { - // Rs has been specified as a function of depth. - OPM_THROW(std::runtime_error, "Cannot initialise: RSVD field not read by EclipseGridParser class."); - } else { - // Default initialisation: constant Rs below contact, saturated above. - for (size_t i = 0; i < rec.size(); ++i) { - const int cell = *(eqlmap.cells(i + 1).begin()); + if (deck.hasField("DISGAS")) { + for (size_t i = 0; i < rec.size(); ++i) { + const int cell = *(eqlmap.cells(i + 1).begin()); + // TODO Check this ... + // The index is here picked as a representative for its equlibrium region, + // but is below used to identify PVT region. + // This assumes that an eq-region has uniform pvt properties, is this always ok? + // For Norne this is trivial, as there is only one active pvt-region (PVTNUM=1 for all cells): + 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" @@ -289,8 +309,40 @@ namespace Opm 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 + 1).begin()); + // TODO Similar as above: eq-region vs pvt-region ... + 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); @@ -311,6 +363,7 @@ namespace Opm typedef EquilReg EqReg; std::vector< std::shared_ptr > rs_func_; + std::vector< std::shared_ptr > rv_func_; PVec pp_; PVec sat_; @@ -335,13 +388,14 @@ namespace Opm const int repcell = *cells.begin(); const RhoCalc calc(props, repcell); - const EqReg eqreg(rec[r], calc, - rs_func_[r], std::make_shared(), + rs_func_[r], rv_func_[r], props.phaseUsage()); - - const PVec press = phasePressures(G, eqreg, cells, grav); + + 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]); @@ -350,8 +404,9 @@ namespace Opm if (props.phaseUsage().phase_used[BlackoilPhases::Liquid] && props.phaseUsage().phase_used[BlackoilPhases::Vapour]) { const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid]; - const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r])); - const Vec rv(cells.size(), 0.0); + 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_); } diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index c075f4ecb..6c9912151 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -580,7 +580,7 @@ namespace Opm phaseSaturations(const Region& reg, const CellRange& cells, const BlackoilPropertiesInterface& props, - const std::vector< std::vector >& phase_pressures) + std::vector< std::vector >& phase_pressures) { const double z0 = reg.datum(); const double zwoc = reg.zwoc (); @@ -621,6 +621,7 @@ namespace Opm 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 @@ -631,8 +632,32 @@ namespace Opm 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; } @@ -659,20 +684,19 @@ namespace Opm std::vector computeRs(const UnstructuredGrid& grid, const CellRangeType& cells, const std::vector oil_pressure, - const Miscibility::RsFunction& rs_func) + 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]); + rs[count] = rs_func(depth, oil_pressure[count], gas_saturation[count]); } return rs; } - - } // namespace Equil From b485d0371199c67ee1d2afd9eb3f5236287fc645 Mon Sep 17 00:00:00 2001 From: osae Date: Fri, 28 Mar 2014 17:18:50 +0100 Subject: [PATCH 49/53] Default equil region should be one not zero ... Otherwise problems when kw EQLNUM is used. --- opm/core/simulator/initStateEquil.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 566f2b00f..cfcb56bae 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -241,8 +241,8 @@ namespace Opm } else { // No explicit equilibration region. - // All cells in region zero. - eqlnum.assign(G.number_of_cells, 0); + // All cells in region one. + eqlnum.assign(G.number_of_cells, 1); } return eqlnum; @@ -384,7 +384,7 @@ namespace Opm r = 0, nr = reg.numRegions(); r < nr; ++r) { - const typename RMap::CellRange cells = reg.cells(r); + const typename RMap::CellRange cells = reg.cells(r+1); const int repcell = *cells.begin(); const RhoCalc calc(props, repcell); From b86fa0af86a898534ba287112feec2dc8fe61825 Mon Sep 17 00:00:00 2001 From: osae Date: Fri, 28 Mar 2014 17:35:43 +0100 Subject: [PATCH 50/53] New parser included. --- examples/compute_initial_state.cpp | 8 +- opm/core/simulator/initStateEquil.hpp | 247 ++++++++++++++++++++- opm/core/simulator/initStateEquil_impl.hpp | 20 ++ 3 files changed, 266 insertions(+), 9 deletions(-) diff --git a/examples/compute_initial_state.cpp b/examples/compute_initial_state.cpp index 0b11ecf1c..3193c9fac 100644 --- a/examples/compute_initial_state.cpp +++ b/examples/compute_initial_state.cpp @@ -29,6 +29,10 @@ #include #include #include + +#include +#include + #include namespace @@ -80,8 +84,10 @@ try 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); + //EclipseGridParser deck(deck_filename); GridManager gm(deck); const UnstructuredGrid& grid = *gm.c_grid(); BlackoilPropertiesFromDeck props(deck, grid, param); diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index cfcb56bae..c4608b804 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -27,6 +27,8 @@ #include #include #include +#include +#include #include #include @@ -63,6 +65,13 @@ namespace Opm 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); /** @@ -227,7 +236,51 @@ namespace Opm } } + 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 @@ -249,6 +302,23 @@ namespace Opm } + inline + std::vector + equilnum(const Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& G ) + { + std::vector eqlnum; + if (newParserDeck->hasKeyword("EQLNUM")) { + eqlnum = newParserDeck->getKeyword("EQLNUM")->getIntData(); + } + else { + // No explicit equilibration region. + // All cells in region one. + eqlnum.assign(G.number_of_cells, 1); + } + + return eqlnum; + } template @@ -278,12 +348,7 @@ namespace Opm rs_func_.reserve(rec.size()); if (deck.hasField("DISGAS")) { for (size_t i = 0; i < rec.size(); ++i) { - const int cell = *(eqlmap.cells(i + 1).begin()); - // TODO Check this ... - // The index is here picked as a representative for its equlibrium region, - // but is below used to identify PVT region. - // This assumes that an eq-region has uniform pvt properties, is this always ok? - // For Norne this is trivial, as there is only one active pvt-region (PVTNUM=1 for all cells): + const int cell = *(eqlmap.cells(i + 1).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 @@ -314,8 +379,7 @@ namespace Opm rv_func_.reserve(rec.size()); if (deck.hasField("VAPOIL")) { for (size_t i = 0; i < rec.size(); ++i) { - const int cell = *(eqlmap.cells(i + 1).begin()); - // TODO Similar as above: eq-region vs pvt-region ... + const int cell = *(eqlmap.cells(i + 1).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 @@ -427,6 +491,173 @@ namespace Opm } }; + + + 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 + 1).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 + 1).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+1); + + 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 diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index 6c9912151..d227d54d1 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -753,6 +753,26 @@ namespace Opm 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. + } From 054e7b42b8fdb4be4b70053c7303ddca7c1165b4 Mon Sep 17 00:00:00 2001 From: osae Date: Mon, 31 Mar 2014 15:32:06 +0200 Subject: [PATCH 51/53] Update tests and provide some eclipse output. --- tests/capillary_overlap.DATA | 99 ++++++++++++++++++++++++++++- tests/equil_liveoil.DATA | 118 +++++++++++++++++++++++++++++++---- tests/test_equil.cpp | 93 ++++++++++++++++++--------- 3 files changed, 268 insertions(+), 42 deletions(-) diff --git a/tests/capillary_overlap.DATA b/tests/capillary_overlap.DATA index 185784eee..b59e42202 100644 --- a/tests/capillary_overlap.DATA +++ b/tests/capillary_overlap.DATA @@ -1,7 +1,77 @@ +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 @@ -22,10 +92,37 @@ SGOF 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 -50 150 50 0.25 45 0.35 +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/equil_liveoil.DATA b/tests/equil_liveoil.DATA index 955a09b7f..e1b417c46 100644 --- a/tests/equil_liveoil.DATA +++ b/tests/equil_liveoil.DATA @@ -1,27 +1,94 @@ +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 / - 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 + 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 @@ -37,10 +104,37 @@ SGOF 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 +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/test_equil.cpp b/tests/test_equil.cpp index b4df1fae8..b1e6f6c0a 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -27,6 +27,9 @@ #include #include +#include +#include + #include #include @@ -339,9 +342,9 @@ BOOST_AUTO_TEST_CASE (DeckAllDead) // solver, and it is unclear if we should check it against // the true answer or something else. const double reltol = 1.0e-3; - BOOST_CHECK_CLOSE(pressures[0][first] , 14955e3 , reltol); - BOOST_CHECK_CLOSE(pressures[0][last ] , 15045e3 , reltol); - BOOST_CHECK_CLOSE(pressures[1][last] , 1.50473e7 , reltol); + BOOST_CHECK_CLOSE(pressures[0][first] , 1.496329839e7 , reltol); + BOOST_CHECK_CLOSE(pressures[0][last ] , 1.50473245e7 , reltol); + BOOST_CHECK_CLOSE(pressures[1][last] , 1.50473245e7 , reltol); } @@ -416,9 +419,9 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary) // solver, and it is unclear if we should check it against // the true answer or something else. const double reltol = 1.0e-6; - BOOST_CHECK_CLOSE(pressures[0][first] , 1.45e7 , reltol); + BOOST_CHECK_CLOSE(pressures[0][first] , 1.469769063e7 , reltol); BOOST_CHECK_CLOSE(pressures[0][last ] , 1.545e7 , reltol); - BOOST_CHECK_CLOSE(pressures[1][last] , 1.5351621345e7 , reltol); + BOOST_CHECK_CLOSE(pressures[1][last] , 1.546e7 , reltol); const auto& sats = comp.saturation(); const std::vector s[3]{ @@ -443,7 +446,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap) Opm::EclipseGridParser deck("capillary_overlap.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 10.0); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -454,9 +457,15 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap) // solver, and it is unclear if we should check it against // the true answer or something else. const double reltol = 1.0e-6; - BOOST_CHECK_CLOSE(pressures[0][first] , 1.45e7 , reltol); - BOOST_CHECK_CLOSE(pressures[0][last ] , 1.545e7 , reltol); - BOOST_CHECK_CLOSE(pressures[1][last] , 1.5351621345e7 , reltol); + const double reltol_ecl = 1.0; + BOOST_CHECK_CLOSE(pressures[0][first], 1.48324e+07, reltol_ecl); // eclipse + BOOST_CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl); + BOOST_CHECK_CLOSE(pressures[1][first], 1.49224e+07, reltol_ecl); + BOOST_CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl); + + BOOST_CHECK_CLOSE(pressures[0][first] , 14832467.14, reltol); // opm + BOOST_CHECK_CLOSE(pressures[0][last ] , 15479883.47, reltol); + BOOST_CHECK_CLOSE(pressures[1][last ] , 15489883.47, reltol); const auto& sats = comp.saturation(); // std::cout << "Saturations:\n"; @@ -466,15 +475,24 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap) // } // std::cout << std::endl; // } - const std::vector s[3]{ - { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.223141818182, 0.532269090909, 0.78471, 0.91526, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.207743333333, 0.08474, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, - { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.776858181818, 0.467730909091, 0.0075466666666, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + + const std::vector s_ecl[3]{// eclipse + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22874042, 0.53397995, 0.78454906, 0.91542006, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20039, 0.08458, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77125955, 0.46602005, 0.015063271, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + + const std::vector s_opm[3]{ // opm + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2289309090909091, 0.53406545454545451, 0.78458, 0.9154, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2002466666666666, 0.0846, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77106909090909093, 0.46593454545454549, 0.015173333333333336, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } }; for (int phase = 0; phase < 3; ++phase) { - BOOST_REQUIRE(sats[phase].size() == s[phase].size()); - for (size_t i = 0; i < s[phase].size(); ++i) { - BOOST_CHECK_CLOSE(sats[phase][i], s[phase][i], reltol); + BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size()); + for (size_t i = 0; i < s_opm[phase].size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl); + BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol); } } } @@ -485,10 +503,12 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil) { Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0); const UnstructuredGrid& grid = *(gm.c_grid()); - Opm::EclipseGridParser deck("equil_liveoil.DATA"); + Opm::ParserPtr parser(new Opm::Parser() ); + Opm::DeckConstPtr deck = parser->parseFile("equil_liveoil.DATA"); + //Opm::EclipseGridParser deck("equil_liveoil.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 10.0); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -499,9 +519,16 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil) // solver, and it is unclear if we should check it against // the true answer or something else. const double reltol = 1.0e-6; - BOOST_CHECK_CLOSE(pressures[0][first], 1.4551302072306179e7, reltol); - BOOST_CHECK_CLOSE(pressures[0][last], 1.5501302072306179e7, reltol); - BOOST_CHECK_CLOSE(pressures[1][last], 1.5538684664272346e7, reltol); + const double reltol_ecl = 1.0; + BOOST_CHECK_CLOSE(pressures[0][first], 1.48324e+07, reltol_ecl); // eclipse + BOOST_CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl); + BOOST_CHECK_CLOSE(pressures[1][first], 1.49224e+07, reltol_ecl); + BOOST_CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl); + + BOOST_CHECK_CLOSE(pressures[0][first], 1.483246714e7, reltol); // opm + BOOST_CHECK_CLOSE(pressures[0][last], 1.547991652e7, reltol); + BOOST_CHECK_CLOSE(pressures[1][first], 1.492246714e7, reltol); + BOOST_CHECK_CLOSE(pressures[1][last], 1.548991652e7, reltol); const auto& sats = comp.saturation(); // std::cout << "Saturations:\n"; @@ -511,16 +538,24 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil) // } // std::cout << std::endl; // } - const std::vector s[3]{ - { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.223141818182, 0.532269090909, 0.78471, 0.91526, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.207743333333, 0.08474, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, - { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.776858181818, 0.467730909091, 0.0075466666666, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + const std::vector s_ecl[3]{ // eclipse + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22898, 0.53422, 0.78470, 0.91531, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20073, 0.08469, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77102, 0.46578, 0.01458, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } }; + + const std::vector s_opm[3]{ // opm + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2291709091, 0.5343054545, 0.78472, 0.91529, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2005866667, 0.08471, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.7708290909, 0.4656945455, 0.01469333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + for (int phase = 0; phase < 3; ++phase) { - BOOST_REQUIRE(sats[phase].size() == s[phase].size()); - for (size_t i = 0; i < s[phase].size(); ++i) { - std::cout << sats[phase][i] << '\n'; - //BOOST_CHECK_CLOSE(sats[phase][i], s[phase][i], reltol); + BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size()); + for (size_t i = 0; i < s_opm[phase].size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol); + BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl); } std::cout << std::endl; } From 3de80989a2d529dabb8766d8917610460d16dfa6 Mon Sep 17 00:00:00 2001 From: osae Date: Mon, 31 Mar 2014 16:16:45 +0200 Subject: [PATCH 52/53] Equil regions now internally indexed 0..(NTEQUL-1) --- opm/core/simulator/initStateEquil.hpp | 34 ++++++++++++++++----------- 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index c4608b804..c1f72b714 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -289,13 +289,16 @@ namespace Opm const UnstructuredGrid& G ) { std::vector eqlnum; - if (deck.hasField("EQLNUM")) { - eqlnum = deck.getIntegerValue("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 one. - eqlnum.assign(G.number_of_cells, 1); + // All cells in region zero. + eqlnum.assign(G.number_of_cells, 0); } return eqlnum; @@ -308,13 +311,16 @@ namespace Opm const UnstructuredGrid& G ) { std::vector eqlnum; - if (newParserDeck->hasKeyword("EQLNUM")) { - eqlnum = newParserDeck->getKeyword("EQLNUM")->getIntData(); + 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 one. - eqlnum.assign(G.number_of_cells, 1); + // All cells in region zero. + eqlnum.assign(G.number_of_cells, 0); } return eqlnum; @@ -348,7 +354,7 @@ namespace Opm rs_func_.reserve(rec.size()); if (deck.hasField("DISGAS")) { for (size_t i = 0; i < rec.size(); ++i) { - const int cell = *(eqlmap.cells(i + 1).begin()); + 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 @@ -379,7 +385,7 @@ namespace Opm rv_func_.reserve(rec.size()); if (deck.hasField("VAPOIL")) { for (size_t i = 0; i < rec.size(); ++i) { - const int cell = *(eqlmap.cells(i + 1).begin()); + 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 @@ -448,7 +454,7 @@ namespace Opm r = 0, nr = reg.numRegions(); r < nr; ++r) { - const typename RMap::CellRange cells = reg.cells(r+1); + const typename RMap::CellRange cells = reg.cells(r); const int repcell = *cells.begin(); const RhoCalc calc(props, repcell); @@ -517,7 +523,7 @@ namespace Opm rs_func_.reserve(rec.size()); if (newParserDeck->hasKeyword("DISGAS")) { for (size_t i = 0; i < rec.size(); ++i) { - const int cell = *(eqlmap.cells(i + 1).begin()); + 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); @@ -547,7 +553,7 @@ namespace Opm rv_func_.reserve(rec.size()); if (newParserDeck->hasKeyword("VAPOIL")) { for (size_t i = 0; i < rec.size(); ++i) { - const int cell = *(eqlmap.cells(i + 1).begin()); + 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); @@ -615,7 +621,7 @@ namespace Opm r = 0, nr = reg.numRegions(); r < nr; ++r) { - const typename RMap::CellRange cells = reg.cells(r+1); + const typename RMap::CellRange cells = reg.cells(r); const int repcell = *cells.begin(); const RhoCalc calc(props, repcell); From a43ae52c18a84f772f65fedac92bd751114f81e8 Mon Sep 17 00:00:00 2001 From: osae Date: Thu, 3 Apr 2014 09:07:00 +0200 Subject: [PATCH 53/53] Some additional tests: live gas, RSVD and RVVD --- tests/capillary.DATA | 9 +- tests/deadfluids.DATA | 9 +- tests/equil_livegas.DATA | 131 ++++++++++++++++++++ tests/equil_rsvd_and_rvvd.DATA | 151 +++++++++++++++++++++++ tests/test_equil.cpp | 219 +++++++++++++++++++++++++++++++-- 5 files changed, 507 insertions(+), 12 deletions(-) create mode 100644 tests/equil_livegas.DATA create mode 100644 tests/equil_rsvd_and_rvvd.DATA diff --git a/tests/capillary.DATA b/tests/capillary.DATA index d8d163fcd..89c938999 100644 --- a/tests/capillary.DATA +++ b/tests/capillary.DATA @@ -2,6 +2,13 @@ WATER OIL GAS +TABDIMS + 1 1 40 20 1 20 / + +EQLDIMS +-- NTEQUL + 1 / + PVDO 100 1.0 1.0 200 0.9 1.0 @@ -27,5 +34,5 @@ DENSITY / EQUIL -50 150 50 0.25 20 0.35 +50 150 50 0.25 20 0.35 1* 1* 0 / diff --git a/tests/deadfluids.DATA b/tests/deadfluids.DATA index ffaebc7ff..89232054b 100644 --- a/tests/deadfluids.DATA +++ b/tests/deadfluids.DATA @@ -2,6 +2,13 @@ WATER OIL GAS +TABDIMS + 1 1 40 20 1 20 / + +EQLDIMS +-- NTEQUL + 1 / + PVDO 100 1.0 1.0 200 0.5 1.0 @@ -27,5 +34,5 @@ DENSITY / EQUIL -5 150 5 0 2 0 +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_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 index b1e6f6c0a..d68470b08 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -329,9 +329,10 @@ BOOST_AUTO_TEST_CASE (DeckAllDead) { std::shared_ptr grid(create_grid_cart3d(1, 1, 10), destroy_grid); - Opm::EclipseGridParser deck("deadfluids.DATA"); + Opm::ParserPtr parser(new Opm::Parser() ); + Opm::DeckConstPtr deck = parser->parseFile("deadfluids.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, *grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, *grid, 10.0); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, *grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells); @@ -354,7 +355,8 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion) // Test setup. Opm::GridManager gm(1, 1, 40, 1.0, 1.0, 2.5); const UnstructuredGrid& grid = *(gm.c_grid()); - Opm::EclipseGridParser deck("capillary.DATA"); + Opm::ParserPtr parser(new Opm::Parser() ); + Opm::DeckConstPtr deck = parser->parseFile("capillary.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); // Test the capillary inversion for oil-water. @@ -405,10 +407,11 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary) { Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0); const UnstructuredGrid& grid = *(gm.c_grid()); - Opm::EclipseGridParser deck("capillary.DATA"); + Opm::ParserPtr parser(new Opm::Parser() ); + Opm::DeckConstPtr deck = parser->parseFile("capillary.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 10.0); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -443,10 +446,11 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap) { Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0); const UnstructuredGrid& grid = *(gm.c_grid()); - Opm::EclipseGridParser deck("capillary_overlap.DATA"); + Opm::ParserPtr parser(new Opm::Parser() ); + Opm::DeckConstPtr deck = parser->parseFile("capillary_overlap.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -505,7 +509,6 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil) const UnstructuredGrid& grid = *(gm.c_grid()); Opm::ParserPtr parser(new Opm::Parser() ); Opm::DeckConstPtr deck = parser->parseFile("equil_liveoil.DATA"); - //Opm::EclipseGridParser deck("equil_liveoil.DATA"); Opm::BlackoilPropertiesFromDeck props(deck, grid, false); Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); @@ -543,13 +546,11 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil) { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20073, 0.08469, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77102, 0.46578, 0.01458, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } }; - const std::vector s_opm[3]{ // opm { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2291709091, 0.5343054545, 0.78472, 0.91529, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2005866667, 0.08471, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.7708290909, 0.4656945455, 0.01469333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } }; - for (int phase = 0; phase < 3; ++phase) { BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size()); for (size_t i = 0; i < s_opm[phase].size(); ++i) { @@ -559,7 +560,205 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil) } std::cout << std::endl; } + + const auto& rs = comp.rs(); + const std::vector rs_opm {74.612335679539058, 74.649052116644228, 74.685786561426298, 74.722539022717172, // opm + 74.759309509353145, 74.796098030174733, 74.8329045940269, 74.869729209758916, + 74.906571886224327, 75.090675116639048, 75.0, 75.0, + 75.0, 75.0, 75.0, 75.0, + 75.0, 75.0, 75.0, 75.0}; + const std::vector rs_ecl {74.612228, 74.648956, 74.685707, 74.722473, // eclipse + 74.759254, 74.796051, 74.832870, 74.875145, + 74.969231, 75.090706, 75.000000, 75.000000, + 75.000000, 75.000000, 75.000000, 75.000000, + 75.000000, 75.000000, 75.000000, 75.000000}; + for (size_t i = 0; i < rs_opm.size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + BOOST_CHECK_CLOSE(rs[i], rs_opm[i], reltol); + BOOST_CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl); + } } + +BOOST_AUTO_TEST_CASE (DeckWithLiveGas) +{ + Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0); + const UnstructuredGrid& grid = *(gm.c_grid()); + Opm::ParserPtr parser(new Opm::Parser() ); + Opm::DeckConstPtr deck = parser->parseFile("equil_livegas.DATA"); + Opm::BlackoilPropertiesFromDeck props(deck, grid, false); + + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); + const auto& pressures = comp.press(); + BOOST_REQUIRE(pressures.size() == 3); + BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-6; + const double reltol_ecl = 1.0; + BOOST_CHECK_CLOSE(pressures[0][first], 1.48215e+07, reltol_ecl); // eclipse + BOOST_CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl); + BOOST_CHECK_CLOSE(pressures[1][first], 1.49115e+07, reltol_ecl); + BOOST_CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl); + + BOOST_CHECK_CLOSE(pressures[0][first], 1.482150311e7, reltol); // opm + BOOST_CHECK_CLOSE(pressures[0][last], 1.547988347e7, reltol); + BOOST_CHECK_CLOSE(pressures[1][first], 1.491150311e7, reltol); + BOOST_CHECK_CLOSE(pressures[1][last], 1.548988347e7, reltol); + + const auto& sats = comp.saturation(); + // std::cout << "Saturations:\n"; + // for (const auto& sat : sats) { + // for (const double s : sat) { + // std::cout << s << ' '; + // } + // std::cout << std::endl; + // } + const std::vector s_ecl[3]{ // eclipse + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.24285614, 0.53869015, 0.78454906, 0.91542006, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18311, 0.08458, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.75714386, 0.46130988, 0.032345835, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + + const std::vector s_opm[3]{ // opm + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.24310545, 0.5388, 0.78458, 0.91540, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18288667, 0.0846, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.75689455, 0.4612, 0.03253333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + for (int phase = 0; phase < 3; ++phase) { + BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size()); + for (size_t i = 0; i < s_opm[phase].size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], 100.*reltol); + BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl); + } + std::cout << std::endl; + } + + const auto& rv = comp.rv(); + const std::vector rv_opm { // opm + 2.4884509e-4, 2.4910378e-4, 2.4936267e-4, 2.4962174e-4, + 2.4988100e-4, 2.5014044e-4, 2.5040008e-4, 2.5065990e-4, + 2.5091992e-4, 2.5118012e-4, 2.5223082e-4, 2.5105e-4, + 2.5105e-4, 2.5105e-4, 2.5105e-4, 2.5105e-4, + 2.5105e-4, 2.5105e-4, 2.5105e-4, 2.5105e-4}; + + const std::vector rv_ecl { // eclipse + 0.24884584E-03, 0.24910446E-03, 0.24936325E-03, 0.24962222E-03, + 0.24988138E-03, 0.25014076E-03, 0.25040031E-03, 0.25066003E-03, + 0.25091995E-03, 0.25118008E-03, 0.25223137E-03, 0.25104999E-03, + 0.25104999E-03, 0.25104999E-03, 0.25104999E-03, 0.25104999E-03, + 0.25104999E-03, 0.25104999E-03, 0.25104999E-03, 0.25104999E-03}; + + for (size_t i = 0; i < rv_opm.size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + BOOST_CHECK_CLOSE(rv[i], rv_opm[i], 100.*reltol); + BOOST_CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl); + } +} + +BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD) +{ + Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0); + const UnstructuredGrid& grid = *(gm.c_grid()); + Opm::ParserPtr parser(new Opm::Parser() ); + Opm::DeckConstPtr deck = parser->parseFile("equil_rsvd_and_rvvd.DATA"); + Opm::BlackoilPropertiesFromDeck props(deck, grid, false); + + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); + const auto& pressures = comp.press(); + BOOST_REQUIRE(pressures.size() == 3); + BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-6; + const double reltol_ecl = 1.0; + BOOST_CHECK_CLOSE(pressures[0][first], 1.48350e+07, reltol_ecl); // eclipse + BOOST_CHECK_CLOSE(pressures[0][last], 1.54794e+07, reltol_ecl); + BOOST_CHECK_CLOSE(pressures[1][first], 1.49250e+07, reltol_ecl); + BOOST_CHECK_CLOSE(pressures[1][last], 1.54894e+07, reltol_ecl); + + BOOST_CHECK_CLOSE(pressures[0][first], 1.483499660e7, reltol); // opm + BOOST_CHECK_CLOSE(pressures[0][last], 1.547924516e7, reltol); + BOOST_CHECK_CLOSE(pressures[1][first], 1.492499660e7, reltol); + BOOST_CHECK_CLOSE(pressures[1][last], 1.548924516e7, reltol); + + const auto& sats = comp.saturation(); + // std::cout << "Saturations:\n"; + // for (const auto& sat : sats) { + // for (const double s : sat) { + // std::cout << s << ' '; + // } + // std::cout << std::endl; + // } + const std::vector s_ecl[3]{ // eclipse + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22206347, 0.52871972, 0.78150368, 0.91819441, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.19656529, 0.081805572, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77793652, 0.47128031, 0.021931054, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + + const std::vector s_opm[3]{ // opm + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22232000, 0.52882909, 0.78153000, 0.91817000, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.19636333, 0.08183000, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77768000, 0.47117091, 0.02210667, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + + for (int phase = 0; phase < 3; ++phase) { + BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size()); + for (size_t i = 0; i < s_opm[phase].size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], 100.*reltol); + BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl); + } + std::cout << std::endl; + } + + const auto& rs = comp.rs(); + const std::vector rs_opm { // opm + 74.624983020822540, 74.659590408801634, 74.694380353364522, 74.729353362649505, + 74.764509945812975, 74.799850613032362, 74.835375875509555, 74.87108624547416, + 74.906982236186707, 75.088917653469309, 52.5, 57.5, + 62.5, 67.5, 72.5, 76.45954840804761, + 76.70621044909619, 76.952877357524045, 77.199549133522638, 77.446225777283587}; + + const std::vector rs_ecl { // eclipse + 74.625114, 74.659706, 74.694481, 74.729439, + 74.764580, 74.799904, 74.835419, 74.875252, + 74.968628, 75.088951, 52.500000, 57.500000, + 62.500000, 67.500000, 72.500000, 76.168388, + 76.349953, 76.531532, 76.713142, 76.894775,}; + + const auto& rv = comp.rv(); + const std::vector rv_opm { // opm + 2.50e-6, 7.50e-6, 1.25e-5, 1.75e-5, + 2.25e-5, 2.75e-5, 3.25e-5, 3.75e-5, + 4.25e-5, 2.51158386e-4, 2.52203372e-4, 5.75e-5, + 6.25e-5, 6.75e-5, 7.25e-5, 7.75e-5, + 8.25e-5, 8.75e-5, 9.25e-5, 9.75e-5}; + + const std::vector rv_ecl { // eclipse + 0.24999999E-05, 0.74999998E-05, 0.12500000E-04, 0.17500000E-04, + 0.22500000E-04, 0.27500000E-04, 0.32500000E-04, 0.37500002E-04, + 0.42500000E-04, 0.25115837E-03, 0.25220393E-03, 0.57500001E-04, + 0.62500003E-04, 0.67499997E-04, 0.72499999E-04, 0.77500001E-04, + 0.82500002E-04, 0.87499997E-04, 0.92499999E-04, 0.97500000E-04}; + + for (size_t i = 0; i < rv_opm.size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + BOOST_CHECK_CLOSE(rs[i], rs_opm[i], 100*reltol); + BOOST_CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl); + BOOST_CHECK_CLOSE(rv[i], rv_opm[i], 100.*reltol); + BOOST_CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl); + } +} + BOOST_AUTO_TEST_SUITE_END()