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).
2014-01-14 13:37:28 -06:00
|
|
|
/*
|
|
|
|
Copyright 2014 SINTEF ICT, Applied Mathematics.
|
|
|
|
|
|
|
|
This file is part of the Open Porous Media project (OPM).
|
|
|
|
|
|
|
|
OPM is free software: you can redistribute it and/or modify
|
|
|
|
it under the terms of the GNU General Public License as published by
|
|
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
|
|
(at your option) any later version.
|
|
|
|
|
|
|
|
OPM is distributed in the hope that it will be useful,
|
|
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
|
GNU General Public License for more details.
|
|
|
|
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
|
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
*/
|
|
|
|
|
|
|
|
#ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
|
|
|
|
#define OPM_INITSTATEEQUIL_HEADER_INCLUDED
|
|
|
|
|
2014-02-24 09:09:04 -06:00
|
|
|
#include <opm/core/simulator/EquilibrationHelpers.hpp>
|
2014-01-23 03:16:49 -06:00
|
|
|
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
|
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).
2014-01-14 13:37:28 -06:00
|
|
|
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
|
|
|
|
#include <opm/core/props/BlackoilPhases.hpp>
|
2014-02-24 08:19:04 -06:00
|
|
|
#include <opm/core/utility/RegionMapping.hpp>
|
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).
2014-01-14 13:37:28 -06:00
|
|
|
#include <opm/core/utility/Units.hpp>
|
|
|
|
|
|
|
|
#include <array>
|
|
|
|
#include <cassert>
|
2014-01-17 13:07:51 -06:00
|
|
|
#include <utility>
|
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).
2014-01-14 13:37:28 -06:00
|
|
|
#include <vector>
|
|
|
|
|
2014-01-19 18:25:33 -06:00
|
|
|
/**
|
|
|
|
* \file
|
|
|
|
* Facilities for an ECLIPSE-style equilibration-based
|
|
|
|
* initialisation scheme (keyword 'EQUIL').
|
|
|
|
*/
|
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).
2014-01-14 13:37:28 -06:00
|
|
|
struct UnstructuredGrid;
|
|
|
|
|
|
|
|
namespace Opm
|
|
|
|
{
|
2014-01-19 18:25:33 -06:00
|
|
|
/**
|
|
|
|
* 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.
|
|
|
|
*/
|
2014-02-24 08:55:14 -06:00
|
|
|
namespace Equil {
|
2014-02-19 06:42:07 -06:00
|
|
|
|
2014-01-19 18:25:33 -06:00
|
|
|
/**
|
|
|
|
* 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
|
2014-01-20 03:33:42 -06:00
|
|
|
* 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.
|
2014-01-19 18:25:33 -06:00
|
|
|
*
|
|
|
|
* \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.
|
|
|
|
*/
|
2014-01-17 10:43:27 -06:00
|
|
|
template <class Region, class CellRange>
|
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).
2014-01-14 13:37:28 -06:00
|
|
|
std::vector< std::vector<double> >
|
|
|
|
phasePressures(const UnstructuredGrid& G,
|
|
|
|
const Region& reg,
|
2014-01-17 10:43:27 -06:00
|
|
|
const CellRange& cells,
|
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).
2014-01-14 13:37:28 -06:00
|
|
|
const double grav = unit::gravity);
|
2014-01-23 03:16:49 -06:00
|
|
|
|
2014-02-19 06:42:07 -06:00
|
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
|
* Compute initial phase saturations by means of equilibration.
|
|
|
|
*
|
|
|
|
* \tparam Region Type of an equilibration region information
|
|
|
|
* base. Typically an instance of the EquilReg
|
|
|
|
* class template.
|
|
|
|
*
|
|
|
|
* \tparam CellRange Type of cell range that demarcates the
|
|
|
|
* cells pertaining to the current
|
|
|
|
* equilibration region. Must implement
|
|
|
|
* methods begin() and end() to bound the range
|
|
|
|
* as well as provide an inner type,
|
|
|
|
* const_iterator, to traverse the range.
|
|
|
|
*
|
|
|
|
* \param[in] reg Current equilibration region.
|
|
|
|
* \param[in] cells Range that spans the cells of the current
|
|
|
|
* equilibration region.
|
|
|
|
* \param[in] props Property object, needed for capillary functions.
|
|
|
|
* \param[in] phase_pressures Phase pressures, one vector for each active phase,
|
|
|
|
* of pressure values in each cell in the current
|
|
|
|
* equilibration region.
|
|
|
|
* \return Phase saturations, one vector for each phase, each containing
|
|
|
|
* one saturation value per cell in the region.
|
|
|
|
*/
|
|
|
|
template <class Region, class CellRange>
|
|
|
|
std::vector< std::vector<double> >
|
|
|
|
phaseSaturations(const Region& reg,
|
|
|
|
const CellRange& cells,
|
|
|
|
const BlackoilPropertiesInterface& props,
|
2014-02-27 02:37:48 -06:00
|
|
|
const std::vector< std::vector<double> >& phase_pressures);
|
2014-02-19 06:42:07 -06:00
|
|
|
|
|
|
|
|
2014-01-23 03:16:49 -06:00
|
|
|
namespace DeckDependent {
|
|
|
|
inline
|
|
|
|
std::vector<EquilRecord>
|
|
|
|
getEquil(const EclipseGridParser& deck)
|
|
|
|
{
|
|
|
|
if (deck.hasField("EQUIL")) {
|
|
|
|
const EQUIL& eql = deck.getEQUIL();
|
|
|
|
|
|
|
|
typedef std::vector<EquilLine>::size_type sz_t;
|
|
|
|
const sz_t nrec = eql.equil.size();
|
|
|
|
|
|
|
|
std::vector<EquilRecord> ret;
|
|
|
|
ret.reserve(nrec);
|
|
|
|
for (sz_t r = 0; r < nrec; ++r) {
|
|
|
|
const EquilLine& rec = eql.equil[r];
|
|
|
|
|
|
|
|
EquilRecord record =
|
|
|
|
{
|
|
|
|
{ rec.datum_depth_ ,
|
|
|
|
rec.datum_depth_pressure_ }
|
|
|
|
,
|
|
|
|
{ rec.water_oil_contact_depth_ ,
|
|
|
|
rec.oil_water_cap_pressure_ }
|
|
|
|
,
|
|
|
|
{ rec.gas_oil_contact_depth_ ,
|
|
|
|
rec.gas_oil_cap_pressure_ }
|
|
|
|
};
|
|
|
|
|
|
|
|
ret.push_back(record);
|
|
|
|
}
|
|
|
|
|
|
|
|
return ret;
|
|
|
|
}
|
|
|
|
else {
|
|
|
|
OPM_THROW(std::domain_error,
|
|
|
|
"Deck does not provide equilibration data.");
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
inline
|
|
|
|
std::vector<int>
|
|
|
|
equilnum(const EclipseGridParser& deck,
|
|
|
|
const UnstructuredGrid& G )
|
|
|
|
{
|
|
|
|
std::vector<int> eqlnum;
|
|
|
|
if (deck.hasField("EQLNUM")) {
|
|
|
|
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 InputDeck>
|
2014-02-27 02:08:39 -06:00
|
|
|
class InitialStateComputer;
|
2014-01-23 03:16:49 -06:00
|
|
|
|
|
|
|
template <>
|
2014-02-27 02:08:39 -06:00
|
|
|
class InitialStateComputer<Opm::EclipseGridParser> {
|
2014-01-23 03:16:49 -06:00
|
|
|
public:
|
2014-02-27 02:08:39 -06:00
|
|
|
InitialStateComputer(const BlackoilPropertiesInterface& props,
|
|
|
|
const EclipseGridParser& deck ,
|
|
|
|
const UnstructuredGrid& G ,
|
|
|
|
const double grav = unit::gravity)
|
2014-01-23 03:16:49 -06:00
|
|
|
: pp_(props.numPhases(),
|
2014-02-21 01:52:25 -06:00
|
|
|
std::vector<double>(G.number_of_cells)),
|
|
|
|
sat_(props.numPhases(),
|
2014-02-27 02:08:39 -06:00
|
|
|
std::vector<double>(G.number_of_cells)),
|
|
|
|
rs_(G.number_of_cells)
|
2014-01-23 03:16:49 -06:00
|
|
|
{
|
2014-02-27 01:31:03 -06:00
|
|
|
// Get the equilibration records.
|
2014-01-23 03:16:49 -06:00
|
|
|
const std::vector<EquilRecord> rec = getEquil(deck);
|
2014-02-27 01:31:03 -06:00
|
|
|
|
|
|
|
// Create (inverse) region mapping.
|
2014-01-23 03:16:49 -06:00
|
|
|
const RegionMapping<> eqlmap(equilnum(deck, G));
|
2014-02-27 01:31:03 -06:00
|
|
|
|
|
|
|
// 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<Miscibility::RsSatAtContact>(props, cell, p_contact));
|
|
|
|
}
|
|
|
|
}
|
|
|
|
} else {
|
|
|
|
for (size_t i = 0; i < rec.size(); ++i) {
|
|
|
|
rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// Compute phase pressures and saturations.
|
2014-02-26 07:49:06 -06:00
|
|
|
calcPressSat(eqlmap, rec, props, G, grav);
|
2014-01-23 03:16:49 -06:00
|
|
|
}
|
|
|
|
|
2014-02-27 02:08:39 -06:00
|
|
|
typedef std::vector<double> Vec;
|
|
|
|
typedef std::vector<Vec> PVec; // One per phase.
|
2014-01-23 03:16:49 -06:00
|
|
|
|
2014-02-27 02:08:39 -06:00
|
|
|
const PVec& press() const { return pp_; }
|
|
|
|
const PVec& saturation() const { return sat_; }
|
2014-01-23 03:16:49 -06:00
|
|
|
|
|
|
|
private:
|
|
|
|
typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
|
|
|
|
typedef EquilReg<RhoCalc> EqReg;
|
|
|
|
|
2014-02-27 01:31:03 -06:00
|
|
|
std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
|
|
|
|
|
2014-02-27 02:08:39 -06:00
|
|
|
PVec pp_;
|
|
|
|
PVec sat_;
|
|
|
|
Vec rs_;
|
2014-01-23 03:16:49 -06:00
|
|
|
|
|
|
|
template <class RMap>
|
|
|
|
void
|
2014-02-26 07:49:06 -06:00
|
|
|
calcPressSat(const RMap& reg ,
|
|
|
|
const std::vector< EquilRecord >& rec ,
|
|
|
|
const Opm::BlackoilPropertiesInterface& props,
|
|
|
|
const UnstructuredGrid& G ,
|
|
|
|
const double grav)
|
2014-02-21 01:52:25 -06:00
|
|
|
{
|
2014-02-24 08:55:14 -06:00
|
|
|
typedef Miscibility::NoMixing NoMix;
|
2014-02-21 01:52:25 -06:00
|
|
|
|
|
|
|
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);
|
|
|
|
|
2014-02-26 07:47:24 -06:00
|
|
|
const EqReg eqreg(rec[r], calc,
|
2014-02-27 01:31:03 -06:00
|
|
|
rs_func_[r], std::make_shared<NoMix>(),
|
2014-02-21 01:52:25 -06:00
|
|
|
props.phaseUsage());
|
|
|
|
|
2014-02-27 02:08:39 -06:00
|
|
|
const PVec press = phasePressures(G, eqreg, cells, grav);
|
|
|
|
const PVec sat = phaseSaturations(eqreg, cells, props, press);
|
|
|
|
const Vec rs(cells.size());// = gasOilRatio();
|
2014-02-21 01:52:25 -06:00
|
|
|
|
2014-02-27 02:31:48 -06:00
|
|
|
const int np = props.numPhases();
|
|
|
|
for (int p = 0; p < np; ++p) {
|
|
|
|
copyFromRegion(press[p], cells, pp_[p]);
|
|
|
|
copyFromRegion(sat[p], cells, sat_[p]);
|
2014-02-21 01:52:25 -06:00
|
|
|
}
|
2014-02-27 02:31:48 -06:00
|
|
|
copyFromRegion(rs, cells, rs_);
|
2014-02-21 01:52:25 -06:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2014-02-27 02:31:48 -06:00
|
|
|
template <class CellRangeType>
|
|
|
|
void copyFromRegion(const Vec& source,
|
|
|
|
const CellRangeType& cells,
|
|
|
|
Vec& destination)
|
|
|
|
{
|
|
|
|
auto s = source.begin();
|
|
|
|
auto c = cells.begin();
|
|
|
|
const auto e = cells.end();
|
|
|
|
for (; c != e; ++c, ++s) {
|
|
|
|
destination[*c] = *s;
|
|
|
|
}
|
|
|
|
}
|
2014-02-21 01:52:25 -06:00
|
|
|
|
2014-01-23 03:16:49 -06:00
|
|
|
};
|
|
|
|
} // namespace DeckDependent
|
2014-02-24 08:55:14 -06:00
|
|
|
} // namespace Equil
|
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).
2014-01-14 13:37:28 -06:00
|
|
|
} // namespace Opm
|
|
|
|
|
|
|
|
#include <opm/core/simulator/initStateEquil_impl.hpp>
|
|
|
|
|
|
|
|
#endif // OPM_INITSTATEEQUIL_HEADER_INCLUDED
|