OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #587 from andlaus/implement_multi-region_pvt_tables

Browse Source 
Implement multi-region PVT properties
This commit is contained in:
Bård Skaflestad 2014-06-06 14:38:36 +02:00
commit 95a4d991f0
39 changed files with 2332 additions and 1901 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
CMakeLists_files.cmake
View File

@ -84,11 +84,11 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/props/pvt/BlackoilPvtProperties.cpp
opm/core/props/pvt/PvtPropertiesBasic.cpp
opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
opm/core/props/pvt/SinglePvtDead.cpp
opm/core/props/pvt/SinglePvtDeadSpline.cpp
opm/core/props/pvt/SinglePvtInterface.cpp
opm/core/props/pvt/SinglePvtLiveGas.cpp
opm/core/props/pvt/SinglePvtLiveOil.cpp
opm/core/props/pvt/PvtDead.cpp
opm/core/props/pvt/PvtDeadSpline.cpp
opm/core/props/pvt/PvtInterface.cpp
opm/core/props/pvt/PvtLiveGas.cpp
opm/core/props/pvt/PvtLiveOil.cpp
opm/core/props/rock/RockBasic.cpp
opm/core/props/rock/RockCompressibility.cpp
opm/core/props/rock/RockFromDeck.cpp
@ -311,12 +311,12 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/props/pvt/BlackoilPvtProperties.hpp
opm/core/props/pvt/PvtPropertiesBasic.hpp
opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp
opm/core/props/pvt/SinglePvtConstCompr.hpp
opm/core/props/pvt/SinglePvtDead.hpp
opm/core/props/pvt/SinglePvtDeadSpline.hpp
opm/core/props/pvt/SinglePvtInterface.hpp
opm/core/props/pvt/SinglePvtLiveGas.hpp
opm/core/props/pvt/SinglePvtLiveOil.hpp
opm/core/props/pvt/PvtConstCompr.hpp
opm/core/props/pvt/PvtDead.hpp
opm/core/props/pvt/PvtDeadSpline.hpp
opm/core/props/pvt/PvtInterface.hpp
opm/core/props/pvt/PvtLiveGas.hpp
opm/core/props/pvt/PvtLiveOil.hpp
opm/core/props/rock/RockBasic.hpp
opm/core/props/rock/RockCompressibility.hpp
opm/core/props/rock/RockFromDeck.hpp

4
opm/core/pressure/CompressibleTpfa.cpp
View File

@ -238,7 +238,7 @@ namespace Opm
const int cell = wells_->well_cells[j];
const double cell_depth = grid_.cell_centroids[dim * cell + dim - 1];
props_.matrix(1, &state.pressure()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0);
props_.density(1, &A[0], &rho[0]);
props_.density(1, &A[0], &cell, &rho[0]);
for (int phase = 0; phase < np; ++phase) {
const double s_phase = state.saturation()[np*cell + phase];
wellperf_wdp_[j] += s_phase*rho[phase]*grav*(cell_depth - ref_depth);
@ -380,7 +380,7 @@ namespace Opm
// Gravity contribution, gravcontrib = rho*(face_z - cell_z) [per phase].
if (grav != 0.0) {
const double depth_diff = face_depth - grid_.cell_centroids[c[j]*dim + dim - 1];
props_.density(1, &cell_A_[np*np*c[j]], &gravcontrib[j][0]);
props_.density(1, &cell_A_[np*np*c[j]], &c[j], &gravcontrib[j][0]);
for (int p = 0; p < np; ++p) {
gravcontrib[j][p] *= depth_diff*grav;
}

4
opm/core/props/BlackoilPropertiesBasic.cpp
View File

@ -157,9 +157,11 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
void BlackoilPropertiesBasic::density(const int n,
const double* A,
const int* /*cells*/,
double* rho) const
{
const int np = numPhases();
@ -177,7 +179,7 @@ namespace Opm
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
const double* BlackoilPropertiesBasic::surfaceDensity() const
const double* BlackoilPropertiesBasic::surfaceDensity(int /*cellIdx*/) const
{
return pvt_.surfaceDensities();
}

10
opm/core/props/BlackoilPropertiesBasic.hpp
View File

@ -59,6 +59,11 @@ namespace Opm
/// \return N, the number of cells.
virtual int numCells() const;
/// Return an array containing the PVT table index for each
/// grid cell
virtual const int* cellPvtRegionIndex() const
{ return NULL; }
/// \return Array of N porosity values.
virtual const double* porosity() const;
@ -114,14 +119,17 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
virtual void density(const int n,
const double* A,
const int* cells,
double* rho) const;
/// Densities of stock components at surface conditions.
/// \param[in] cellIdx The index of the cell for which the surface density ought to be calculated
/// \return Array of P density values.
virtual const double* surfaceDensity() const;
virtual const double* surfaceDensity(int cellIdx = 0) const;
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.

37
opm/core/props/BlackoilPropertiesFromDeck.cpp
View File

@ -96,14 +96,20 @@ namespace Opm
void BlackoilPropertiesFromDeck::viscosity(const int n,
const double* p,
const double* z,
const int* /*cells*/,
const int* cells,
double* mu,
double* dmudp) const
{
if (dmudp) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::viscosity() -- derivatives of viscosity not yet implemented.");
} else {
pvt_.mu(n, p, z, mu);
const int *cellPvtTableIdx = cellPvtRegionIndex();
assert(cellPvtTableIdx != 0);
std::vector<int> pvtTableIdx(n);
for (int i = 0; i < n; ++ i)
pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
pvt_.mu(n, &pvtTableIdx[0], p, z, mu);
}
}
@ -120,21 +126,27 @@ namespace Opm
void BlackoilPropertiesFromDeck::matrix(const int n,
const double* p,
const double* z,
const int* /*cells*/,
const int* cells,
double* A,
double* dAdp) const
{
const int np = numPhases();
const int *cellPvtTableIdx = cellPvtRegionIndex();
std::vector<int> pvtTableIdx(n);
for (int i = 0; i < n; ++ i)
pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
B_.resize(n*np);
R_.resize(n*np);
if (dAdp) {
dB_.resize(n*np);
dR_.resize(n*np);
pvt_.dBdp(n, p, z, &B_[0], &dB_[0]);
pvt_.dRdp(n, p, z, &R_[0], &dR_[0]);
pvt_.dBdp(n, &pvtTableIdx[0], p, z, &B_[0], &dB_[0]);
pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]);
} else {
pvt_.B(n, p, z, &B_[0]);
pvt_.R(n, p, z, &R_[0]);
pvt_.B(n, &pvtTableIdx[0], p, z, &B_[0]);
pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]);
}
const int* phase_pos = pvt_.phasePosition();
bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] &&
@ -207,15 +219,19 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
void BlackoilPropertiesFromDeck::density(const int n,
const double* A,
const int* cells,
double* rho) const
{
const int np = numPhases();
const double* sdens = pvt_.surfaceDensities();
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int cellIdx = cells?cells[i]:i;
int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx);
const double* sdens = pvt_.surfaceDensities(pvtRegionIdx);
for (int phase = 0; phase < np; ++phase) {
rho[np*i + phase] = 0.0;
for (int comp = 0; comp < np; ++comp) {
@ -227,9 +243,10 @@ namespace Opm
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
const double* BlackoilPropertiesFromDeck::surfaceDensity() const
const double* BlackoilPropertiesFromDeck::surfaceDensity(int cellIdx) const
{
return pvt_.surfaceDensities();
int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx);
return pvt_.surfaceDensities(pvtRegionIdx);
}
/// \param[in] n Number of data points.

17
opm/core/props/BlackoilPropertiesFromDeck.hpp
View File

@ -97,6 +97,11 @@ namespace Opm
/// \return N, the number of cells.
virtual int numCells() const;
/// Return an array containing the PVT table index for each
/// grid cell
virtual const int* cellPvtRegionIndex() const
{ return &cellPvtRegionIdx_[0]; }
/// \return Array of N porosity values.
virtual const double* porosity() const;
@ -152,14 +157,16 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
virtual void density(const int n,
const double* A,
const int* cells,
double* rho) const;
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
virtual const double* surfaceDensity() const;
virtual const double* surfaceDensity(int cellIdx = 0) const;
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
@ -206,6 +213,13 @@ namespace Opm
double* smax) const;
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
template<class CentroidIterator>
void init(Opm::DeckConstPtr deck,
int number_of_cells,
@ -224,6 +238,7 @@ namespace Opm
const parameter::ParameterGroup& param,
bool init_rock);
RockFromDeck rock_;
std::vector<int> cellPvtRegionIdx_;
BlackoilPvtProperties pvt_;
std::unique_ptr<SaturationPropsInterface> satprops_;
mutable std::vector<double> B_;

17
opm/core/props/BlackoilPropertiesFromDeck_impl.hpp
View File

@ -23,8 +23,14 @@ namespace Opm
const parameter::ParameterGroup& param,
bool init_rock)
{
init(deck, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
param, init_rock);
init(deck,
number_of_cells,
global_cell,
cart_dims,
begin_cell_centroids,
dimension,
param,
init_rock);
}
template<class CentroidIterator>
@ -36,6 +42,10 @@ namespace Opm
int dimension,
bool init_rock)
{
// retrieve the cell specific PVT table index from the deck
// and using the grid...
extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell);
if (init_rock){
rock_.init(deck, number_of_cells, global_cell, cart_dims);
}
@ -62,6 +72,9 @@ namespace Opm
const parameter::ParameterGroup& param,
bool init_rock)
{
// retrieve the cell specific PVT table index from the deck
// and using the grid...
extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell);
if(init_rock){
rock_.init(deck, number_of_cells, global_cell, cart_dims);

8
opm/core/props/BlackoilPropertiesInterface.hpp
View File

@ -47,6 +47,10 @@ namespace Opm
/// \return N, the number of cells.
virtual int numCells() const = 0;
/// Return an array containing the PVT table index for each
/// grid cell
virtual const int* cellPvtRegionIndex() const = 0;
/// \return Array of N porosity values.
virtual const double* porosity() const = 0;
@ -102,14 +106,16 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
virtual void density(const int n,
const double* A,
const int* cells,
double* rho) const = 0;
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
virtual const double* surfaceDensity() const = 0;
virtual const double* surfaceDensity(int regionIdx = 0) const = 0;
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.

120
opm/core/props/pvt/BlackoilPvtProperties.cpp
View File

@ -21,11 +21,11 @@
#include "config.h"
#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
#include <opm/core/props/pvt/SinglePvtDead.hpp>
#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
#include <opm/core/props/pvt/SinglePvtConstCompr.hpp>
#include <opm/core/props/pvt/PvtDead.hpp>
#include <opm/core/props/pvt/PvtDeadSpline.hpp>
#include <opm/core/props/pvt/PvtLiveOil.hpp>
#include <opm/core/props/pvt/PvtLiveGas.hpp>
#include <opm/core/props/pvt/PvtConstCompr.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
@ -42,93 +42,97 @@ namespace Opm
{
}
void BlackoilPvtProperties::init(Opm::DeckConstPtr deck, int samples)
void BlackoilPvtProperties::init(Opm::DeckConstPtr deck,
int numSamples)
{
// If we need multiple regions, this class and the SinglePvt* classes must change.
region_number_ = 0;
phase_usage_ = phaseUsageFromDeck(deck);
// Surface densities. Accounting for different orders in eclipse and our code.
if (deck->hasKeyword("DENSITY")) {
Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY");
int numRegions = densityKeyword->size();
densities_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
if (phase_usage_.phase_used[Liquid]) {
densities_[phase_usage_.phase_pos[Liquid]]
= densityKeyword->getRecord(region_number_)->getItem("OIL")->getSIDouble(0);
densities_[regionIdx][phase_usage_.phase_pos[Liquid]]
= densityKeyword->getRecord(regionIdx)->getItem("OIL")->getSIDouble(0);
}
if (phase_usage_.phase_used[Aqua]) {
densities_[phase_usage_.phase_pos[Aqua]]
= densityKeyword->getRecord(region_number_)->getItem("WATER")->getSIDouble(0);
densities_[regionIdx][phase_usage_.phase_pos[Aqua]]
= densityKeyword->getRecord(regionIdx)->getItem("WATER")->getSIDouble(0);
}
if (phase_usage_.phase_used[Vapour]) {
densities_[phase_usage_.phase_pos[Vapour]]
= densityKeyword->getRecord(region_number_)->getItem("GAS")->getSIDouble(0);
densities_[regionIdx][phase_usage_.phase_pos[Vapour]]
= densityKeyword->getRecord(regionIdx)->getItem("GAS")->getSIDouble(0);
}
} else {
OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");
}
// Set the properties.
// Resize the property objects container
props_.resize(phase_usage_.num_phases);
// Water PVT
if (phase_usage_.phase_used[Aqua]) {
if (deck->hasKeyword("PVTW")) {
Opm::PvtwTable pvtwTable(deck->getKeyword("PVTW"));
// if water is used, we require the presence of the "PVTW"
// keyword for now...
std::shared_ptr<PvtConstCompr> pvtw(new PvtConstCompr);
pvtw->initFromWater(deck->getKeyword("PVTW"));
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable));
} else {
// Eclipse 100 default.
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
}
props_[phase_usage_.phase_pos[Aqua]] = pvtw;
}
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
// for oil, we support the "PVDO", "PVTO" and "PVCDO"
// keywords...
if (deck->hasKeyword("PVDO")) {
Opm::PvdoTable pvdoTable(deck->getKeyword("PVDO"), region_number_);
if (samples > 0) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
Opm::DeckKeywordConstPtr pvdoKeyword = deck->getKeyword("PVDO");
if (numSamples > 0) {
auto splinePvt = std::shared_ptr<PvtDeadSpline>(new PvtDeadSpline);
splinePvt->initFromOil(pvdoKeyword, numSamples);
props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
} else {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable));
auto deadPvt = std::shared_ptr<PvtDead>(new PvtDead);
deadPvt->initFromOil(pvdoKeyword);
props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
}
} else if (deck->hasKeyword("PVTO")) {
Opm::PvtoTable pvtoTable(deck->getKeyword("PVTO"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable));
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO")));
} else if (deck->hasKeyword("PVCDO")) {
Opm::PvcdoTable pvcdoTable(deck->getKeyword("PVCDO"));
std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
pvcdo->initFromOil(deck->getKeyword("PVCDO"));
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable));
props_[phase_usage_.phase_pos[Liquid]] = pvcdo;
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n");
}
}
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
// gas can be specified using the "PVDG" or "PVTG" keywords...
if (deck->hasKeyword("PVDG")) {
Opm::PvdgTable pvdgTable(deck->getKeyword("PVDG"), region_number_);
if (samples > 0) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
Opm::DeckKeywordConstPtr pvdgKeyword = deck->getKeyword("PVDG");
if (numSamples > 0) {
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromGas(pvdgKeyword, numSamples);
props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
} else {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable));
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromGas(pvdgKeyword);
props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
}
} else if (deck->hasKeyword("PVTG")) {
Opm::PvtgTable pvtgTable(deck->getKeyword("PVTG"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable));
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG")));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
}
}
// Must inform pvt property objects of phase structure.
for (int i = 0; i < phase_usage_.num_phases; ++i) {
props_[i]->setPhaseConfiguration(phase_usage_.num_phases, phase_usage_.phase_pos);
}
}
const double* BlackoilPvtProperties::surfaceDensities() const
const double* BlackoilPvtProperties::surfaceDensities(int regionIdx) const
{
return densities_;
return &densities_[regionIdx][0];
}
@ -154,13 +158,14 @@ namespace Opm
void BlackoilPvtProperties::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const
{
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->mu(n, p, z, &data1_[0]);
props_[phase]->mu(n, pvtTableIdx, p, z, &data1_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[phase_usage_.num_phases*i + phase] = data1_[i];
@ -169,13 +174,14 @@ namespace Opm
}
void BlackoilPvtProperties::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B) const
{
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->B(n, p, z, &data1_[0]);
props_[phase]->B(n, pvtTableIdx, p, z, &data1_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i];
@ -184,6 +190,7 @@ namespace Opm
}
void BlackoilPvtProperties::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@ -192,7 +199,7 @@ namespace Opm
data1_.resize(n);
data2_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->dBdp(n, p, z, &data1_[0], &data2_[0]);
props_[phase]->dBdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i];
@ -203,13 +210,14 @@ namespace Opm
void BlackoilPvtProperties::R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const
{
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->R(n, p, z, &data1_[0]);
props_[phase]->R(n, pvtTableIdx, p, z, &data1_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[phase_usage_.num_phases*i + phase] = data1_[i];
@ -218,6 +226,7 @@ namespace Opm
}
void BlackoilPvtProperties::dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
@ -226,7 +235,7 @@ namespace Opm
data1_.resize(n);
data2_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->dRdp(n, p, z, &data1_[0], &data2_[0]);
props_[phase]->dRdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[phase_usage_.num_phases*i + phase] = data1_[i];
@ -234,5 +243,4 @@ namespace Opm
}
}
}
} // namespace Opm

35
opm/core/props/pvt/BlackoilPvtProperties.hpp
View File

@ -20,26 +20,30 @@
#ifndef OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
#define OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
#include <string>
#include <memory>
#include <array>
namespace Opm
{
/// Class collecting the pvt properties for all active phases.
/// For all the methods, the following apply: p and z
/// are expected to be of size n and n*num_phases, respectively.
/// Output arrays shall be of size n*num_phases, and must be valid
/// For all the methods, the following apply:
/// - p and z are expected to be of size n and n*num_phases, respectively.
/// - pvtTableIdx specifies the PVT table to be used for each data
/// point and is thus expected to be an array of size n
/// - Output arrays shall be of size n*num_phases, and must be valid
/// before calling the method.
/// NOTE: The difference between this interface and the one defined
/// by SinglePvtInterface is that this collects all phases' properties,
/// by PvtInterface is that this collects all phases' properties,
/// and therefore the output arrays are of size n*num_phases as opposed
/// to size n in SinglePvtInterface.
/// to size n in PvtInterface.
class BlackoilPvtProperties : public BlackoilPhases
{
public:
@ -47,8 +51,10 @@ namespace Opm
BlackoilPvtProperties();
/// Initialize from deck.
///
/// \param deck An input deck from the opm-parser module.
void init(Opm::DeckConstPtr deck, int samples);
void init(Opm::DeckConstPtr deck,
int samples);
/// \return Object describing the active phases.
PhaseUsage phaseUsage() const;
@ -68,22 +74,25 @@ namespace Opm
/// Densities of stock components at surface conditions.
/// \return Array of size numPhases().
const double* surfaceDensities() const;
const double* surfaceDensities(int regionIdx = 0) const;
/// Viscosity as a function of p and z.
void mu(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const;
/// Formation volume factor as a function of p and z.
void B(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
void dBdp(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@ -91,12 +100,14 @@ namespace Opm
/// Solution factor as a function of p and z.
void R(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
void dRdp(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_R,
@ -109,11 +120,11 @@ namespace Opm
PhaseUsage phase_usage_;
int region_number_;
// The PVT properties. We need to store one object per PVT
// region per active fluid phase.
std::vector<std::shared_ptr<PvtInterface> > props_;
std::vector<std::array<double, MaxNumPhases> > densities_;
std::vector<std::shared_ptr<SinglePvtInterface> > props_;
double densities_[MaxNumPhases];
mutable std::vector<double> data1_;
mutable std::vector<double> data2_;
};

307
opm/core/props/pvt/PvtConstCompr.hpp Normal file
View File

@ -0,0 +1,307 @@
/*
Copyright 2010, 2011, 2012 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_PVTCONSTCOMPR_HEADER_INCLUDED
#define OPM_PVTCONSTCOMPR_HEADER_INCLUDED
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
#include <opm/parser/eclipse/Utility/PvcdoTable.hpp>
#include <vector>
#include <algorithm>
namespace Opm
{
/// Class for constant compressible phases (PVTW or PVCDO).
/// The PVT properties can either be given as a function of
/// pressure (p) and surface volume (z) or pressure (p) and gas
/// resolution factor (r). Also, since this class supports
/// multiple PVT regions, the concrete table to be used for each
/// data point needs to be specified via the pvtTableIdx argument
/// of the respective method. For all the virtual methods, the
/// following apply: pvtTableIdx, p, r and z are expected to be of
/// size n, size n, size n and n*num_phases, respectively. Output
/// arrays shall be of size n, and must be valid before calling
/// the method.
class PvtConstCompr : public PvtInterface
{
public:
PvtConstCompr()
{}
void initFromWater(Opm::DeckKeywordConstPtr pvtwKeyword)
{
int numRegions = pvtwKeyword->size();
ref_press_.resize(numRegions);
ref_B_.resize(numRegions);
comp_.resize(numRegions);
viscosity_.resize(numRegions);
visc_comp_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
Opm::DeckRecordConstPtr pvtwRecord = pvtwKeyword->getRecord(regionIdx);
ref_press_[regionIdx] = pvtwRecord->getItem("P_REF")->getSIDouble(0);
ref_B_[regionIdx] = pvtwRecord->getItem("WATER_VOL_FACTOR")->getSIDouble(0);
comp_[regionIdx] = pvtwRecord->getItem("WATER_COMPRESSIBILITY")->getSIDouble(0);
viscosity_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSITY")->getSIDouble(0);
visc_comp_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSIBILITY")->getSIDouble(0);
}
}
void initFromOil(Opm::DeckKeywordConstPtr pvcdoKeyword)
{
int numRegions = pvcdoKeyword->size();
ref_press_.resize(numRegions);
ref_B_.resize(numRegions);
comp_.resize(numRegions);
viscosity_.resize(numRegions);
visc_comp_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
Opm::DeckRecordConstPtr pvcdoRecord = pvcdoKeyword->getRecord(regionIdx);
ref_press_[regionIdx] = pvcdoRecord->getItem("P_REF")->getSIDouble(0);
ref_B_[regionIdx] = pvcdoRecord->getItem("OIL_VOL_FACTOR")->getSIDouble(0);
comp_[regionIdx] = pvcdoRecord->getItem("OIL_COMPRESSIBILITY")->getSIDouble(0);
viscosity_[regionIdx] = pvcdoRecord->getItem("OIL_VISCOSITY")->getSIDouble(0);
visc_comp_[regionIdx] = pvcdoRecord->getItem("OIL_VISCOSIBILITY")->getSIDouble(0);
}
}
/*!
* \brief Create a PVT object with a given viscosity that
* assumes all fluid phases to be incompressible.
*/
explicit PvtConstCompr(double visc)
: ref_press_(1, 0.0),
ref_B_(1, 1.0),
comp_(1, 0.0),
viscosity_(1, visc),
visc_comp_(1, 0.0)
{
}
virtual ~PvtConstCompr()
{
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*z*/,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
output_mu[i] = viscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
}
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*r*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
output_mu[i] = viscosity_[tableIdx]/d;
output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx];
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
output_mu[i] = viscosity_[tableIdx]/d;
output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx];
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*z*/,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
output_B[i] = ref_B_[tableIdx]/(1.0 + x + 0.5*x*x);
}
}
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
output_B[i] = ref_B_[tableIdx]/d;
output_dBdp[i] = (-ref_B_[tableIdx]/(d*d))*(1 + x) * comp_[tableIdx];
}
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*r*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
// b = 1/B = d/ref_B_B;
output_b[i] = d/ref_B_[tableIdx];
output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx];
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
// b = 1/B = d/ref_B_[tableIdx]B;
output_b[i] = d/ref_B_[tableIdx];
output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx];
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
virtual void rsSat(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
virtual void rvSat(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
{
std::fill(output_rvSat, output_rvSat + n, 0.0);
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
virtual void R(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
{
std::fill(output_R, output_R + n, 0.0);
}
virtual void dRdp(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R,
double* output_dRdp) const
{
std::fill(output_R, output_R + n, 0.0);
std::fill(output_dRdp, output_dRdp + n, 0.0);
}
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
// The PVT properties. We need to store one value per PVT
// region.
std::vector<double> ref_press_;
std::vector<double> ref_B_;
std::vector<double> comp_;
std::vector<double> viscosity_;
std::vector<double> visc_comp_;
};
}
#endif // OPM_PVTCONSTCOMPR_HEADER_INCLUDED

117
opm/core/props/pvt/SinglePvtDead.cpp → opm/core/props/pvt/PvtDead.cpp
View File

@ -19,7 +19,7 @@
#include "config.h"
#include <opm/core/props/pvt/SinglePvtDead.hpp>
#include <opm/core/props/pvt/PvtDead.hpp>
#include <algorithm>
// Extra includes for debug dumping of tables.
@ -34,8 +34,16 @@ namespace Opm
// Member functions
//-------------------------------------------------------------------------
/// Constructor
SinglePvtDead::SinglePvtDead(const Opm::PvdoTable& pvdoTable)
void PvtDead::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword)
{
int numRegions = Opm::PvdoTable::numTables(pvdoKeyword);
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx);
// Copy data
const std::vector<double>& press = pvdoTable.getPressureColumn();
const std::vector<double>& b = pvdoTable.getFormationFactorColumn();
@ -46,20 +54,23 @@ namespace Opm
for (int i = 0; i < sz; ++i) {
bInv[i] = 1.0 / b[i];
}
b_ = NonuniformTableLinear<double>(press, bInv);
viscosity_ = NonuniformTableLinear<double>(press, visc);
// Dumping the created tables.
// static int count = 0;
// std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
// os.precision(15);
// os << "1/B\n\n" << one_over_B_
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
b_[regionIdx] = NonuniformTableLinear<double>(press, bInv);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
}
}
/// Constructor
SinglePvtDead::SinglePvtDead(const Opm::PvdgTable& pvdgTable)
void PvtDead::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword)
{
int numRegions = Opm::PvdgTable::numTables(pvdgKeyword);
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx);
// Copy data
const std::vector<double>& press = pvdgTable.getPressureColumn();
const std::vector<double>& b = pvdgTable.getFormationFactorColumn();
@ -70,36 +81,33 @@ namespace Opm
for (int i = 0; i < sz; ++i) {
bInv[i] = 1.0 / b[i];
}
b_ = NonuniformTableLinear<double>(press, bInv);
viscosity_ = NonuniformTableLinear<double>(press, visc);
// Dumping the created tables.
// static int count = 0;
// std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
// os.precision(15);
// os << "1/B\n\n" << one_over_B_
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
b_[regionIdx] = NonuniformTableLinear<double>(press, bInv);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
}
}
// Destructor
SinglePvtDead::~SinglePvtDead()
PvtDead::~PvtDead()
{
}
void SinglePvtDead::mu(const int n,
void PvtDead::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*z*/,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = viscosity_(p[i]);
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
}
}
void SinglePvtDead::mu(const int n,
void PvtDead::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
double* output_mu,
@ -108,14 +116,16 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = viscosity_(p[i]);
output_dmudp[i] = viscosity_.derivative(p[i]);
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void SinglePvtDead::mu(const int n,
void PvtDead::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
@ -125,14 +135,16 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = viscosity_(p[i]);
output_dmudp[i] = viscosity_.derivative(p[i]);
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void SinglePvtDead::B(const int n,
void PvtDead::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*z*/,
double* output_B) const
@ -140,25 +152,29 @@ namespace Opm
// #pragma omp parallel for
// B = 1/b
for (int i = 0; i < n; ++i) {
output_B[i] = 1.0/b_(p[i]);
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_B[i] = 1.0/b_[regionIdx](p[i]);
}
}
void SinglePvtDead::dBdp(const int n,
void PvtDead::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
B(n, p, 0, output_B);
B(n, pvtTableIdx, p, 0, output_B);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
double Bg = output_B[i];
output_dBdp[i] = -Bg*Bg*b_.derivative(p[i]);
output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]);
}
}
void SinglePvtDead::b(const int n,
void PvtDead::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
double* output_b,
@ -168,15 +184,18 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_b[i] = b_(p[i]);
output_dbdp[i] = b_.derivative(p[i]);
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void SinglePvtDead::b(const int n,
void PvtDead::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
@ -187,15 +206,18 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_b[i] = b_(p[i]);
output_dbdp[i] = b_.derivative(p[i]);
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void SinglePvtDead::rsSat(const int n,
void PvtDead::rsSat(const int n,
const int* pvtTableIdx,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
@ -204,7 +226,8 @@ namespace Opm
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
void SinglePvtDead::rvSat(const int n,
void PvtDead::rvSat(const int n,
const int* pvtTableIdx,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
@ -213,7 +236,8 @@ namespace Opm
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
void SinglePvtDead::R(const int n,
void PvtDead::R(const int n,
const int* pvtTableIdx,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
@ -221,7 +245,8 @@ namespace Opm
std::fill(output_R, output_R + n, 0.0);
}
void SinglePvtDead::dRdp(const int n,
void PvtDead::dRdp(const int n,
const int* pvtTableIdx,
const double* /*p*/,
const double* /*z*/,
double* output_R,

60
opm/core/props/pvt/SinglePvtDead.hpp → opm/core/props/pvt/PvtDead.hpp
View File

@ -17,11 +17,11 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SINGLEPVTDEAD_HEADER_INCLUDED
#define OPM_SINGLEPVTDEAD_HEADER_INCLUDED
#ifndef OPM_PVTDEAD_HEADER_INCLUDED
#define OPM_PVTDEAD_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
@ -33,21 +33,28 @@ namespace Opm
{
/// Class for immiscible dead oil and dry gas.
/// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
/// or pressure (p) and gas resolution factor (r).
/// For all the virtual methods, the following apply: p, r and z
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
/// calling the method.
class SinglePvtDead : public SinglePvtInterface
/// The PVT properties can either be given as a function of
/// pressure (p) and surface volume (z) or pressure (p) and gas
/// resolution factor (r). Also, since this class supports
/// multiple PVT regions, the concrete table to be used for each
/// data point needs to be specified via the pvtTableIdx argument
/// of the respective method. For all the virtual methods, the
/// following apply: pvtTableIdx, p, r and z are expected to be of
/// size n, size n, size n and n*num_phases, respectively. Output
/// arrays shall be of size n, and must be valid before calling
/// the method.
class PvtDead : public PvtInterface
{
public:
SinglePvtDead(const Opm::PvdoTable& pvdoTable);
SinglePvtDead(const Opm::PvdgTable& pvdgTable);
virtual ~SinglePvtDead();
PvtDead() {};
void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword);
void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword);
virtual ~PvtDead();
/// Viscosity as a function of p and z.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const;
@ -55,6 +62,7 @@ namespace Opm
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_mu,
@ -64,6 +72,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@ -73,12 +82,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@ -87,6 +98,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_b,
@ -96,6 +108,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@ -106,35 +119,46 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
private:
// PVT properties of dry gas or dead oil
NonuniformTableLinear<double> b_;
NonuniformTableLinear<double> viscosity_;
};
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
// PVT properties of dry gas or dead oil. We need to store one
// table per PVT region.
std::vector<NonuniformTableLinear<double> > b_;
std::vector<NonuniformTableLinear<double> > viscosity_;
};
}
#endif // OPM_SINGLEPVTDEAD_HEADER_INCLUDED
#endif // OPM_PVTDEAD_HEADER_INCLUDED

260
opm/core/props/pvt/PvtDeadSpline.cpp Normal file
View File

@ -0,0 +1,260 @@
/*
Copyright 2010, 2011, 2012 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/>.
*/
#include "config.h"
#include <opm/core/props/pvt/PvtDeadSpline.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/parser/eclipse/Utility/SingleRecordTable.hpp>
#include <algorithm>
// Extra includes for debug dumping of tables.
// #include <boost/lexical_cast.hpp>
// #include <string>
// #include <fstream>
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
PvtDeadSpline::PvtDeadSpline()
{}
void PvtDeadSpline::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
int numSamples)
{
int numRegions = Opm::PvdoTable::numTables(pvdoKeyword);
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx);
int numRows = pvdoTable.numRows();
// Copy data
const std::vector<double> &press = pvdoTable.getPressureColumn();
const std::vector<double> &B = pvdoTable.getFormationFactorColumn();
const std::vector<double> &visc = pvdoTable.getViscosityColumn();
std::vector<double> B_inv(numRows);
for (int i = 0; i < numRows; ++i) {
B_inv[i] = 1.0 / B[i];
}
buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]);
buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]);
}
}
void PvtDeadSpline::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
int numSamples)
{
int numRegions = Opm::PvdgTable::numTables(pvdgKeyword);
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx);
int numRows = pvdgTable.numRows();
// Copy data
const std::vector<double> &press = pvdgTable.getPressureColumn();
const std::vector<double> &B = pvdgTable.getFormationFactorColumn();
const std::vector<double> &visc = pvdgTable.getViscosityColumn();
std::vector<double> B_inv(numRows);
for (int i = 0; i < numRows; ++i) {
B_inv[i] = 1.0 / B[i];
}
buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]);
buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]);
}
}
// Destructor
PvtDeadSpline::~PvtDeadSpline()
{
}
void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*z*/,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
}
}
void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void PvtDeadSpline::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*z*/,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_B[i] = 1.0/b_[regionIdx](p[i]);
}
}
void PvtDeadSpline::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
B(n, pvtTableIdx, p, 0, output_B);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
double Bg = output_B[i];
output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]);
}
}
void PvtDeadSpline::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void PvtDeadSpline::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void PvtDeadSpline::rsSat(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
void PvtDeadSpline::rvSat(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
{
std::fill(output_rvSat, output_rvSat + n, 0.0);
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
void PvtDeadSpline::R(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
{
std::fill(output_R, output_R + n, 0.0);
}
void PvtDeadSpline::dRdp(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R,
double* output_dRdp) const
{
std::fill(output_R, output_R + n, 0.0);
std::fill(output_dRdp, output_dRdp + n, 0.0);
}
}

46
opm/core/props/pvt/SinglePvtDeadSpline.hpp → opm/core/props/pvt/PvtDeadSpline.hpp
View File

@ -17,10 +17,10 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
#define OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
#ifndef OPM_PVTDEADSPLINE_HEADER_INCLUDED
#define OPM_PVTDEADSPLINE_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
@ -38,15 +38,21 @@ namespace Opm
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
/// calling the method.
class SinglePvtDeadSpline : public SinglePvtInterface
class PvtDeadSpline : public PvtInterface
{
public:
SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples);
SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples);
virtual ~SinglePvtDeadSpline();
PvtDeadSpline();
void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
int numSamples);
void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
int numSamples);
virtual ~PvtDeadSpline();
/// Viscosity as a function of p and z.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const;
@ -54,6 +60,7 @@ namespace Opm
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_mu,
@ -63,6 +70,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@ -72,12 +80,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@ -86,6 +96,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_b,
@ -95,6 +106,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@ -104,35 +116,47 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
private:
// PVT properties of dry gas or dead oil
UniformTableLinear<double> b_;
UniformTableLinear<double> viscosity_;
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
// PVT properties of dry gas or dead oil. We need to store one
// table per PVT region.
std::vector<UniformTableLinear<double> > b_;
std::vector<UniformTableLinear<double> > viscosity_;
};
}
#endif // OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
#endif // OPM_PVTDEADSPLINE_HEADER_INCLUDED

79
opm/core/props/pvt/PvtInterface.cpp Normal file
View File

@ -0,0 +1,79 @@
/*
Copyright 2010, 2011, 2012 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/>.
*/
#include "config.h"
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <cassert>
namespace Opm
{
PvtInterface::PvtInterface()
: num_phases_(MaxNumPhases)
{
for (int i = 0; i < MaxNumPhases; ++i) {
phase_pos_[i] = i;
}
}
PvtInterface::~PvtInterface()
{
}
void PvtInterface::setPhaseConfiguration(const int num_phases, const int* phase_pos)
{
num_phases_ = num_phases;
for (int i = 0; i < MaxNumPhases; ++i) {
phase_pos_[i] = phase_pos[i];
}
}
void extractPvtTableIndex(std::vector<int> &pvtTableIdx,
Opm::DeckConstPtr deck,
size_t numCompressed,
const int *compressedToCartesianCellIdx)
{
if (deck->hasKeyword("PVTNUM")) {
// we have a "real" multi-pvt deck.
// get the PVT number from the deck
Opm::DeckKeywordConstPtr pvtnumKeyword = deck->getKeyword("PVTNUM");
const std::vector<int> &pvtnumData = pvtnumKeyword->getIntData();
// convert this into an array of compressed cells. since
// Eclipse uses Fortran-style indices which start at 1
// instead of 0, we subtract 1.
pvtTableIdx.resize(numCompressed);
for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) {
size_t cartesianCellIdx = compressedToCartesianCellIdx?compressedToCartesianCellIdx[cellIdx]:cellIdx;
assert(cartesianCellIdx < pvtnumData.size());
pvtTableIdx[cellIdx] = pvtnumData[cartesianCellIdx] - 1;
}
}
else {
// create the pvtIdxArray: all cells use the one and only
// PVT table...
pvtTableIdx.resize(numCompressed);
std::fill(pvtTableIdx.begin(), pvtTableIdx.end(), 0);
}
}
} // namespace Opm

54
opm/core/props/pvt/SinglePvtInterface.hpp → opm/core/props/pvt/PvtInterface.hpp
View File

@ -17,22 +17,22 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED
#define OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED
#ifndef OPM_PVTINTERFACE_HEADER_INCLUDED
#define OPM_PVTINTERFACE_HEADER_INCLUDED
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
namespace Opm
{
class SinglePvtInterface : public BlackoilPhases
class PvtInterface : public BlackoilPhases
{
public:
SinglePvtInterface();
PvtInterface();
virtual ~SinglePvtInterface();
virtual ~PvtInterface();
/// \param[in] num_phases The number of active phases.
/// \param[in] phase_pos Array of BlackpoilPhases::MaxNumPhases
@ -44,13 +44,19 @@ namespace Opm
/// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
/// or pressure (p) and gas resolution factor (r).
/// For all the virtual methods, the following apply: p, r and z
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
/// For all the virtual methods, the following apply:
/// - pvtRegionIdx is an array of size n and represents the
/// index of the PVT table which should be used to calculate
/// the output. NULL can also be passed and is interpreted
/// such that the first table should be used for the output
/// - p, r and z are expected to be of size n, size n and
/// n*num_phases, respectively.
/// - Output arrays shall be of size n, and must be valid before
/// calling the method.
/// Viscosity as a function of p and z.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_mu) const = 0;
@ -58,6 +64,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
double* output_mu,
@ -67,6 +74,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@ -76,12 +84,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B) const = 0;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B,
@ -90,6 +100,7 @@ namespace Opm
/// The inverse of the volume factor b = 1 / B as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
double* output_b,
@ -99,6 +110,7 @@ namespace Opm
/// The inverse of the volume factor b = 1 / B as a function of p and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@ -108,12 +120,14 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const = 0;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const = 0;
@ -121,12 +135,14 @@ namespace Opm
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const = 0;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
@ -138,7 +154,25 @@ namespace Opm
int phase_pos_[MaxNumPhases];
};
/*!
* \brief Helper function to create an array containing the (C-Style)
* PVT table index for each compressed cell from an Eclipse deck.
*
* This function assumes that the degrees of freedom where PVT
* properties need to be calculated are grid cells. The main point
* of this function is to avoid code duplication because the
* Eclipse deck only contains Fortran-style PVT table indices
* which start at 1 instead of 0 and -- more relevantly -- it uses
* logically cartesian cell indices to specify the table index of
* a cell while the classes which use the PvtInterface
* implementations usually use compressed cells.
*/
void extractPvtTableIndex(std::vector<int>& pvtTableIdx,
Opm::DeckConstPtr deck,
size_t numCompressed,
const int* compressedToCartesianIdx);
} // namespace Opm
#endif // OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED
#endif // OPM_PVTINTERFACE_HEADER_INCLUDED

526
opm/core/props/pvt/PvtLiveGas.cpp Normal file
View File

@ -0,0 +1,526 @@
//===========================================================================
//
// File: MiscibilityLiveGas.cpp
//
// Created: Wed Feb 10 09:21:53 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
/*
Copyright 2010 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/>.
*/
#include "config.h"
#include <opm/core/props/pvt/PvtLiveGas.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <algorithm>
#include <iostream>
namespace Opm
{
using Opm::linearInterpolation;
using Opm::linearInterpolationDerivative;
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
PvtLiveGas::PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword)
{
int numTables = Opm::PvtgTable::numTables(pvtgKeyword);
saturated_gas_table_.resize(numTables);
undersat_gas_tables_.resize(numTables);
for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) {
Opm::PvtgTable pvtgTable(pvtgKeyword, pvtTableIdx);
// GAS, PVTG
saturated_gas_table_[pvtTableIdx].resize(4);
saturated_gas_table_[pvtTableIdx][0] = pvtgTable.getOuterTable()->getPressureColumn();
saturated_gas_table_[pvtTableIdx][1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn();
saturated_gas_table_[pvtTableIdx][2] = pvtgTable.getOuterTable()->getGasViscosityColumn();
saturated_gas_table_[pvtTableIdx][3] = pvtgTable.getOuterTable()->getOilSolubilityColumn();
int sz = pvtgTable.getOuterTable()->numRows();
undersat_gas_tables_[pvtTableIdx].resize(sz);
for (int i=0; i<sz; ++i) {
const auto &undersatTable = *pvtgTable.getInnerTable(i);
undersat_gas_tables_[pvtTableIdx][i].resize(3);
undersat_gas_tables_[pvtTableIdx][i][0] = undersatTable.getOilSolubilityColumn();
undersat_gas_tables_[pvtTableIdx][i][1] = undersatTable.getGasFormationFactorColumn();
undersat_gas_tables_[pvtTableIdx][i][2] = undersatTable.getGasViscosityColumn();
}
// Bg -> 1/Bg
for (int i=0; i<sz; ++i) {
saturated_gas_table_[pvtTableIdx][1][i] = 1.0/saturated_gas_table_[pvtTableIdx][1][i];
for (size_t j=0; j<undersat_gas_tables_[pvtTableIdx][i][1].size(); ++j) {
undersat_gas_tables_[pvtTableIdx][i][1][j] = 1.0/undersat_gas_tables_[pvtTableIdx][i][1][j];
}
}
}
}
// Destructor
PvtLiveGas::~PvtLiveGas()
{
}
void PvtLiveGas::mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 2, false);
}
}
/// Viscosity and its derivatives as a function of p and r.
void PvtLiveGas::mu(const int /*n*/,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*r*/,
double* /*output_mu*/,
double* /*output_dmudp*/,
double* /*output_dmudr*/) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
/// Viscosity and its derivatives as a function of p and r.
void PvtLiveGas::mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
int tableIdx = getTableIndex_(pvtRegionIdx, i);
output_mu[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 0);
output_dmudp[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 1);
output_dmudr[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 2);
}
}
/// Formation volume factor as a function of p and z.
void PvtLiveGas::B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[i] = evalB(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i));
}
}
/// Formation volume factor and p-derivative as functions of p and z.
void PvtLiveGas::dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B,
double* output_dBdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalBDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_B[i], output_dBdp[i]);
}
}
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
void PvtLiveGas::b(const int /*n*/,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*r*/,
double* /*output_b*/,
double* /*output_dbdp*/,
double* /*output_dbdr*/) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
void PvtLiveGas::b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
int tableIdx = getTableIndex_(pvtRegionIdx, i);
output_b[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 0);
output_dbdp[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 1);
output_dbdr[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 2);
}
}
/// Gas resolution and its derivatives at bublepoint as a function of p.
void PvtLiveGas::rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const
{
for (int i = 0; i < n; ++i) {
int pvtTableIdx = getTableIndex_(pvtRegionIdx, i);
output_rvSat[i] = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][3],p[i]);
output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][3],p[i]);
}
}
void PvtLiveGas::rsSat(const int n,
const int* pvtRegionIdx,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
/// Solution factor as a function of p and z.
void PvtLiveGas::R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[i] = evalR(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i));
}
}
/// Solution factor and p-derivative as functions of p and z.
void PvtLiveGas::dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalRDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_R[i], output_dRdp[i]);
}
}
// ---- Private methods ----
double PvtLiveGas::evalB(const double press, const double* surfvol, int pvtTableIdx) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
// To handle no-gas case.
return 1.0;
}
return 1.0/miscible_gas(press, surfvol, pvtTableIdx, 1, false);
}
void PvtLiveGas::evalBDeriv(const double press, const double* surfvol, int pvtTableIdx,
double& Bval, double& dBdpval) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
// To handle no-gas case.
Bval = 1.0;
dBdpval = 0.0;
return;
}
Bval = evalB(press, surfvol, pvtTableIdx);
dBdpval = -Bval*Bval*miscible_gas(press, surfvol, pvtTableIdx, 1, true);
}
double PvtLiveGas::evalR(const double press, const double* surfvol, int pvtTableIdx) const
{
if (surfvol[phase_pos_[Liquid]] == 0.0) {
// To handle no-gas case.
return 0.0;
}
double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][3], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) {
// Saturated case
return satR;
} else {
// Undersaturated case
return maxR;
}
}
void PvtLiveGas::evalRDeriv(const double press, const double* surfvol, int pvtTableIdx,
double& Rval, double& dRdpval) const
{
if (surfvol[phase_pos_[Liquid]] == 0.0) {
// To handle no-gas case.
Rval = 0.0;
dRdpval = 0.0;
return;
}
double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][3], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) {
// Saturated case
Rval = satR;
dRdpval = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][3],
press);
} else {
// Undersaturated case
Rval = maxR;
dRdpval = 0.0;
}
}
double PvtLiveGas::miscible_gas(const double press,
const double* surfvol,
const int pvtTableIdx,
const int item,
const bool deriv) const
{
const std::vector<std::vector<double> > &saturatedGasTable =
saturated_gas_table_[pvtTableIdx];
const std::vector<std::vector<std::vector<double> > > &undersatGasTables =
undersat_gas_tables_[pvtTableIdx];
int section;
double Rval = linearInterpolation(saturatedGasTable[0],
saturatedGasTable[3], press,
section);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (deriv) {
if (Rval < maxR ) { // Saturated case
return linearInterpolationDerivative(saturatedGasTable[0],
saturatedGasTable[item],
press);
} else { // Undersaturated case
int is = section;
if (undersatGasTables[is][0].size() < 2) {
double val = (saturatedGasTable[item][is+1]
- saturatedGasTable[item][is]) /
(saturatedGasTable[0][is+1] -
saturatedGasTable[0][is]);
return val;
}
double val1 =
linearInterpolation(undersatGasTables[is][0],
undersatGasTables[is][item],
maxR);
double val2 =
linearInterpolation(undersatGasTables[is+1][0],
undersatGasTables[is+1][item],
maxR);
double val = (val2 - val1)/
(saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
return val;
}
} else {
if (Rval < maxR ) { // Saturated case
return linearInterpolation(saturatedGasTable[0],
saturatedGasTable[item],
press);
} else { // Undersaturated case
int is = section;
// Extrapolate from first table section
if (is == 0 && press < saturatedGasTable[0][0]) {
return linearInterpolation(undersatGasTables[0][0],
undersatGasTables[0][item],
maxR);
}
// Extrapolate from last table section
int ltp = saturatedGasTable[0].size() - 1;
if (is+1 == ltp && press > saturatedGasTable[0][ltp]) {
return linearInterpolation(undersatGasTables[ltp][0],
undersatGasTables[ltp][item],
maxR);
}
// Interpolate between table sections
double w = (press - saturatedGasTable[0][is]) /
(saturatedGasTable[0][is+1] -
saturatedGasTable[0][is]);
if (undersatGasTables[is][0].size() < 2) {
double val = saturatedGasTable[item][is] +
w*(saturatedGasTable[item][is+1] -
saturatedGasTable[item][is]);
return val;
}
double val1 =
linearInterpolation(undersatGasTables[is][0],
undersatGasTables[is][item],
maxR);
double val2 =
linearInterpolation(undersatGasTables[is+1][0],
undersatGasTables[is+1][item],
maxR);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
double PvtLiveGas::miscible_gas(const double press,
const double r,
const PhasePresence& cond,
const int pvtTableIdx,
const int item,
const int deriv) const
{
const std::vector<std::vector<double> > &saturatedGasTable =
saturated_gas_table_[pvtTableIdx];
const std::vector<std::vector<std::vector<double> > > &undersatGasTables =
undersat_gas_tables_[pvtTableIdx];
const bool isSat = cond.hasFreeOil();
// Derivative w.r.t p
if (deriv == 1) {
if (isSat) { // Saturated case
return linearInterpolationDerivative(saturatedGasTable[0],
saturatedGasTable[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturatedGasTable[0], press);
if (undersatGasTables[is][0].size() < 2) {
double val = (saturatedGasTable[item][is+1]
- saturatedGasTable[item][is]) /
(saturatedGasTable[0][is+1] -
saturatedGasTable[0][is]);
return val;
}
double val1 =
linearInterpolation(undersatGasTables[is][0],
undersatGasTables[is][item],
r);
double val2 =
linearInterpolation(undersatGasTables[is+1][0],
undersatGasTables[is+1][item],
r);
double val = (val2 - val1)/
(saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
return val;
}
} else if (deriv == 2){
if (isSat) {
return 0;
} else {
int is = tableIndex(saturatedGasTable[0], press);
double w = (press - saturatedGasTable[0][is]) /
(saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
assert(undersatGasTables[is][0].size() >= 2);
assert(undersatGasTables[is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersatGasTables[is][0],
undersatGasTables[is][item],
r);
double val2 =
linearInterpolationDerivative(undersatGasTables[is+1][0],
undersatGasTables[is+1][item],
r);
double val = val1 + w * (val2 - val1);
return val;
}
} else {
if (isSat) { // Saturated case
return linearInterpolation(saturatedGasTable[0],
saturatedGasTable[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturatedGasTable[0], press);
// Extrapolate from first table section
if (is == 0 && press < saturatedGasTable[0][0]) {
return linearInterpolation(undersatGasTables[0][0],
undersatGasTables[0][item],
r);
}
// Extrapolate from last table section
//int ltp = saturatedGasTable[0].size() - 1;
//if (is+1 == ltp && press > saturatedGasTable[0][ltp]) {
// return linearInterpolation(undersatGasTables[ltp][0],
// undersatGasTables[ltp][item],
// r);
//}
// Interpolate between table sections
double w = (press - saturatedGasTable[0][is]) /
(saturatedGasTable[0][is+1] -
saturatedGasTable[0][is]);
if (undersatGasTables[is][0].size() < 2) {
double val = saturatedGasTable[item][is] +
w*(saturatedGasTable[item][is+1] -
saturatedGasTable[item][is]);
return val;
}
double val1 =
linearInterpolation(undersatGasTables[is][0],
undersatGasTables[is][item],
r);
double val2 =
linearInterpolation(undersatGasTables[is+1][0],
undersatGasTables[is+1][item],
r);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

50
opm/core/props/pvt/SinglePvtLiveGas.hpp → opm/core/props/pvt/PvtLiveGas.hpp
View File

@ -17,10 +17,10 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
#define OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
#ifndef OPM_PVTLIVEGAS_HEADER_INCLUDED
#define OPM_PVTLIVEGAS_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/parser/eclipse/Utility/PvtgTable.hpp>
@ -35,14 +35,15 @@ namespace Opm
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
/// calling the method.
class SinglePvtLiveGas : public SinglePvtInterface
class PvtLiveGas : public PvtInterface
{
public:
SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable);
virtual ~SinglePvtLiveGas();
PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword);
virtual ~PvtLiveGas();
/// Viscosity as a function of p and z.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_mu) const;
@ -50,6 +51,7 @@ namespace Opm
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
double* output_mu,
@ -59,6 +61,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@ -68,12 +71,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B,
@ -82,6 +87,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
double* output_b,
@ -91,6 +97,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@ -102,52 +109,65 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
protected:
double evalB(double press, const double* surfvol) const;
void evalBDeriv(double press, const double* surfvol, double& B, double& dBdp) const;
double evalR(double press, const double* surfvol) const;
void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const;
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
double evalB(double press, const double* surfvol, int pvtTableIdx) const;
void evalBDeriv(double press, const double* surfvol, int pvtTableIdx, double& B, double& dBdp) const;
double evalR(double press, const double* surfvol, int pvtTableIdx) const;
void evalRDeriv(double press, const double* surfvol, int pvtTableIdx, double& R, double& dRdp) const;
// item: 1=>B 2=>mu;
double miscible_gas(const double press,
const double* surfvol,
const int pvtTableIdx,
const int item,
const bool deriv = false) const;
double miscible_gas(const double press,
const double r,
const PhasePresence& cond,
const int pvtTableIdx,
const int item,
const int deriv = 0) const;
// PVT properties of wet gas (with vaporised oil)
std::vector<std::vector<double> > saturated_gas_table_;
std::vector<std::vector<std::vector<double> > > undersat_gas_tables_;
// PVT properties of wet gas (with vaporised oil). We need to
// store one table per PVT region.
std::vector< std::vector<std::vector<double> > > saturated_gas_table_;
std::vector< std::vector<std::vector<std::vector<double> > > > undersat_gas_tables_;
};
}
#endif // OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
#endif // OPM_PVTLIVEGAS_HEADER_INCLUDED

591
opm/core/props/pvt/PvtLiveOil.cpp Normal file
View File

@ -0,0 +1,591 @@
/*
Copyright 2010, 2011, 2012 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/>.
*/
#include "config.h"
#include <opm/core/props/pvt/PvtLiveOil.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <algorithm>
namespace Opm
{
using Opm::linearInterpolation;
using Opm::linearInterpolationDerivative;
using Opm::tableIndex;
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
PvtLiveOil::PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword)
{
int numTables = Opm::PvtoTable::numTables(pvtoKeyword);
saturated_oil_table_.resize(numTables);
undersat_oil_tables_.resize(numTables);
for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) {
Opm::PvtoTable pvtoTable(pvtoKeyword, pvtTableIdx);
const auto saturatedPvto = pvtoTable.getOuterTable();
// OIL, PVTO
saturated_oil_table_[pvtTableIdx].resize(4);
const int sz = saturatedPvto->numRows();
for (int k=0; k<4; ++k) {
saturated_oil_table_[pvtTableIdx][k].resize(sz);
}
for (int i=0; i<sz; ++i) {
saturated_oil_table_[pvtTableIdx][0][i] = saturatedPvto->getPressureColumn()[i]; // p
saturated_oil_table_[pvtTableIdx][1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo
saturated_oil_table_[pvtTableIdx][2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o
saturated_oil_table_[pvtTableIdx][3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs
}
undersat_oil_tables_[pvtTableIdx].resize(sz);
for (int i=0; i<sz; ++i) {
const auto undersaturatedPvto = pvtoTable.getInnerTable(i);
undersat_oil_tables_[pvtTableIdx][i].resize(3);
int tsize = undersaturatedPvto->numRows();
undersat_oil_tables_[pvtTableIdx][i][0].resize(tsize);
undersat_oil_tables_[pvtTableIdx][i][1].resize(tsize);
undersat_oil_tables_[pvtTableIdx][i][2].resize(tsize);
for (int j=0; j<tsize; ++j) {
undersat_oil_tables_[pvtTableIdx][i][0][j] = undersaturatedPvto->getPressureColumn()[j]; // p
undersat_oil_tables_[pvtTableIdx][i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j]; // 1/Bo
undersat_oil_tables_[pvtTableIdx][i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j]; // mu_o
}
}
// Complete undersaturated tables by extrapolating from existing data
// as is done in Eclipse and Mrst
int iNext = -1;
for (int i=0; i<sz; ++i) {
// Skip records already containing undersaturated data
if (undersat_oil_tables_[pvtTableIdx][i][0].size() > 1) {
continue;
}
// Look ahead for next record containing undersaturated data
if (iNext < i) {
iNext = i+1;
while (iNext<sz && undersat_oil_tables_[pvtTableIdx][iNext][0].size() < 2) {
++iNext;
}
if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
}
// Add undersaturated data to current record while maintaining compressibility and viscosibility
typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
for (sz_t j=1; j<undersat_oil_tables_[pvtTableIdx][iNext][0].size(); ++j) {
double diffPressure = undersat_oil_tables_[pvtTableIdx][iNext][0][j]-undersat_oil_tables_[pvtTableIdx][iNext][0][j-1];
double pressure = undersat_oil_tables_[pvtTableIdx][i][0].back()+diffPressure;
undersat_oil_tables_[pvtTableIdx][i][0].push_back(pressure);
double compr = (1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j]-1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j-1])
/ (0.5*(1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j]+1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j-1]));
double B = (1.0/undersat_oil_tables_[pvtTableIdx][i][1].back())*(1.0+0.5*compr)/(1.0-0.5*compr);
undersat_oil_tables_[pvtTableIdx][i][1].push_back(1.0/B);
double visc = (undersat_oil_tables_[pvtTableIdx][iNext][2][j]-undersat_oil_tables_[pvtTableIdx][iNext][2][j-1])
/ (0.5*(undersat_oil_tables_[pvtTableIdx][iNext][2][j]+undersat_oil_tables_[pvtTableIdx][iNext][2][j-1]));
double mu = (undersat_oil_tables_[pvtTableIdx][i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc);
undersat_oil_tables_[pvtTableIdx][i][2].push_back(mu);
}
}
}
}
/// Destructor.
PvtLiveOil::~PvtLiveOil()
{
}
/// Viscosity as a function of p and z.
void PvtLiveOil::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = miscible_oil(p[i], z + num_phases_*i, tableIdx, 2, false);
}
}
/// Viscosity and its derivatives as a function of p and r.
void PvtLiveOil::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = miscible_oil(p[i], r[i], tableIdx, 2, 0);
output_dmudp[i] = miscible_oil(p[i], r[i], tableIdx, 2, 1);
output_dmudr[i] = miscible_oil(p[i], r[i], tableIdx, 2, 2);
}
}
/// Viscosity and its derivatives as a function of p and r.
void PvtLiveOil::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
const PhasePresence& cnd = cond[i];
output_mu[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 0);
output_dmudp[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 1);
output_dmudr[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 2);
}
}
/// Formation volume factor as a function of p and z.
void PvtLiveOil::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_B[i] = evalB(tableIdx, p[i], z + num_phases_*i);
}
}
/// Formation volume factor and p-derivative as functions of p and z.
void PvtLiveOil::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B,
double* output_dBdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
evalBDeriv(tableIdx, p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
}
}
void PvtLiveOil::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = miscible_oil(tableIdx, p[i], r[i], 1, 0);
output_dbdp[i] = miscible_oil(tableIdx, p[i], r[i], 1, 1);
output_dbdr[i] = miscible_oil(tableIdx, p[i], r[i], 1, 2);
}
}
void PvtLiveOil::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 0);
output_dbdp[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 1);
output_dbdr[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 2);
}
}
void PvtLiveOil::rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const
{
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_rsSat[i] = linearInterpolation(saturated_oil_table_[tableIdx][0],
saturated_oil_table_[tableIdx][3],p[i]);
output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[tableIdx][0],
saturated_oil_table_[tableIdx][3],p[i]);
}
}
void PvtLiveOil::rvSat(const int n,
const int* pvtTableIdx,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
{
std::fill(output_rvSat, output_rvSat + n, 0.0);
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
/// Solution factor as a function of p and z.
void PvtLiveOil::R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_R[i] = evalR(tableIdx, p[i], z + num_phases_*i);
}
}
/// Solution factor and p-derivative as functions of p and z.
void PvtLiveOil::dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
evalRDeriv(tableIdx, p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
}
}
// ---- Private methods ----
double PvtLiveOil::evalB(size_t pvtTableIdx,
double press, const double* surfvol) const
{
// if (surfvol[phase_pos_[Liquid]] == 0.0) return 1.0; // To handle no-oil case.
return 1.0/miscible_oil(press, surfvol, pvtTableIdx, 1, false);
}
void PvtLiveOil::evalBDeriv(size_t pvtTableIdx,
const double press, const double* surfvol,
double& Bval, double& dBdpval) const
{
Bval = evalB(pvtTableIdx, press, surfvol);
dBdpval = -Bval*Bval*miscible_oil(press, surfvol, pvtTableIdx, 1, true);
}
double PvtLiveOil::evalR(size_t pvtTableIdx,
double press, const double* surfvol) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
return 0.0;
}
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][3], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) { // Saturated case
return Rval;
} else {
return maxR; // Undersaturated case
}
}
void PvtLiveOil::evalRDeriv(size_t pvtTableIdx,
const double press, const double* surfvol,
double& Rval, double& dRdpval) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
Rval = 0.0;
dRdpval = 0.0;
return;
}
Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][3], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) {
// Saturated case
dRdpval = linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][3],
press);
} else {
// Undersaturated case
Rval = maxR;
dRdpval = 0.0;
}
}
double PvtLiveOil::miscible_oil(const double press,
const double* surfvol,
const int pvtTableIdx,
const int item,
const bool deriv) const
{
int section;
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][3],
press, section);
double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (deriv) {
if (Rval < maxR ) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], maxR);
double w = (maxR - saturated_oil_table_[pvtTableIdx][3][is]) /
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
} else {
if (Rval < maxR ) { // Saturated case
return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], maxR);
double w = (maxR - saturated_oil_table_[pvtTableIdx][3][is]) /
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
double PvtLiveOil::miscible_oil(const double press,
const double r,
const int pvtTableIdx,
const int item,
const int deriv) const
{
int section;
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][3],
press, section);
// derivative with respect to frist component (pressure)
if (deriv == 1) {
if (Rval <= r ) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
// derivative with respect to second component (r)
} else if (deriv == 2) {
if (Rval <= r ) { // Saturated case
return 0;
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][3][is+1]-saturated_oil_table_[pvtTableIdx][3][is]);
return val;
}
} else {
if (Rval <= r ) { // Saturated case
return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
double PvtLiveOil::miscible_oil(const double press,
const double r,
const PhasePresence& cond,
const int pvtTableIdx,
const int item,
const int deriv) const
{
const bool isSat = cond.hasFreeGas();
// derivative with respect to frist component (pressure)
if (deriv == 1) {
if (isSat) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
// derivative with respect to second component (r)
} else if (deriv == 2) {
if (isSat) { // Saturated case
return 0;
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][3][is+1]-saturated_oil_table_[pvtTableIdx][3][is]);
return val;
}
} else {
if (isSat) { // Saturated case
return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

51
opm/core/props/pvt/SinglePvtLiveOil.hpp → opm/core/props/pvt/PvtLiveOil.hpp
View File

@ -17,11 +17,12 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
#define OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
#ifndef OPM_PVTLIVEOIL_HEADER_INCLUDED
#define OPM_PVTLIVEOIL_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
#include <vector>
@ -35,14 +36,15 @@ namespace Opm
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
/// calling the method.
class SinglePvtLiveOil : public SinglePvtInterface
class PvtLiveOil : public PvtInterface
{
public:
SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable);
virtual ~SinglePvtLiveOil();
PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword);
virtual ~PvtLiveOil();
/// Viscosity as a function of p and z.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const;
@ -50,6 +52,7 @@ namespace Opm
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_mu,
@ -59,6 +62,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@ -68,12 +72,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@ -82,6 +88,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_b,
@ -91,6 +98,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@ -100,58 +108,73 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
private:
double evalB(double press, const double* surfvol) const;
void evalBDeriv(double press, const double* surfvol, double& B, double& dBdp) const;
double evalR(double press, const double* surfvol) const;
void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const;
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
double evalB(size_t pvtTableIdx, double press, const double* surfvol) const;
void evalBDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& B, double& dBdp) const;
double evalR(size_t pvtTableIdx, double press, const double* surfvol) const;
void evalRDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& R, double& dRdp) const;
// item: 1=>1/B 2=>mu;
double miscible_oil(const double press,
const double* surfvol,
const int pvtTableIdx,
const int item,
const bool deriv = false) const;
double miscible_oil(const double press,
const double r,
const int pvtTableIdx,
const int item,
const int deriv = 0) const;
double miscible_oil(const double press,
const double r,
const PhasePresence& cond,
const int pvtTableIdx,
const int item,
const int deriv = 0) const;
// PVT properties of live oil (with dissolved gas)
std::vector<std::vector<double> > saturated_oil_table_;
std::vector<std::vector<std::vector<double> > > undersat_oil_tables_;
// PVT properties of live oil (with dissolved gas). We need to
// store one table per PVT region.
std::vector<std::vector<std::vector<double> > > saturated_oil_table_;
std::vector<std::vector<std::vector<std::vector<double> > > > undersat_oil_tables_;
};
}
#endif // OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
#endif // OPM_PVTLIVEOIL_HEADER_INCLUDED

2
opm/core/props/pvt/PvtPropertiesBasic.hpp
View File

@ -95,6 +95,8 @@ namespace Opm
double* output_dRdp) const;
private:
// The PVT properties. We need to store one value per PVT
// region.
std::vector<double> density_;
std::vector<double> viscosity_;
std::vector<double> formation_volume_factor_;

2
opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
View File

@ -35,7 +35,7 @@ namespace Opm
void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck)
{
// If we need multiple regions, this class and the SinglePvt* classes must change.
// So far, this class only supports a single PVT region. TODO?
int region_number = 0;
PhaseUsage phase_usage = phaseUsageFromDeck(deck);

291
opm/core/props/pvt/SinglePvtConstCompr.hpp
View File

@ -1,291 +0,0 @@
/*
Copyright 2010, 2011, 2012 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_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
#define OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
#include <opm/parser/eclipse/Utility/PvcdoTable.hpp>
#include <vector>
#include <algorithm>
namespace Opm
{
/// Class for constant compressible phases (PVTW or PVCDO).
/// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
/// or pressure (p) and gas resolution factor (r).
/// For all the virtual methods, the following apply: p, r and z
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
/// calling the method.
class SinglePvtConstCompr : public SinglePvtInterface
{
public:
SinglePvtConstCompr(const Opm::PvtwTable &pvtwTable)
{
if (pvtwTable.numRows() != 1)
OPM_THROW(std::runtime_error,
"The table specified by the PVTW keyword is required"
"to exhibit exactly one row.");
ref_press_ = pvtwTable.getPressureColumn()[0];
ref_B_ = pvtwTable.getFormationFactorColumn()[0];
comp_ = pvtwTable.getCompressibilityColumn()[0];
viscosity_ = pvtwTable.getViscosityColumn()[0];
visc_comp_ = pvtwTable.getViscosibilityColumn()[0];
}
SinglePvtConstCompr(const Opm::PvcdoTable &pvcdoTable)
{
if (pvcdoTable.numRows() != 1)
OPM_THROW(std::runtime_error,
"The table specified by the PVCDO keyword is required"
"to exhibit exactly one row.");
ref_press_ = pvcdoTable.getPressureColumn()[0];
ref_B_ = pvcdoTable.getFormationFactorColumn()[0];
comp_ = pvcdoTable.getCompressibilityColumn()[0];
viscosity_ = pvcdoTable.getViscosityColumn()[0];
visc_comp_ = pvcdoTable.getViscosibilityColumn()[0];
}
SinglePvtConstCompr(double visc)
: ref_press_(0.0),
ref_B_(1.0),
comp_(0.0),
viscosity_(visc),
visc_comp_(0.0)
{
}
virtual ~SinglePvtConstCompr()
{
}
virtual void mu(const int n,
const double* p,
const double* /*z*/,
double* output_mu) const
{
if (visc_comp_) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
double x = -visc_comp_*(p[i] - ref_press_);
output_mu[i] = viscosity_/(1.0 + x + 0.5*x*x);
}
} else {
std::fill(output_mu, output_mu + n, viscosity_);
}
}
virtual void mu(const int n,
const double* p,
const double* /*r*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
if (visc_comp_) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
double x = -visc_comp_*(p[i] - ref_press_);
double d = (1.0 + x + 0.5*x*x);
output_mu[i] = viscosity_/d;
output_dmudp[i] = (viscosity_/(d*d))*(1+x) * visc_comp_;
}
} else {
std::fill(output_mu, output_mu + n, viscosity_);
std::fill(output_dmudp, output_dmudp + n, 0.0);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
virtual void mu(const int n,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
if (visc_comp_) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
double x = -visc_comp_*(p[i] - ref_press_);
double d = (1.0 + x + 0.5*x*x);
output_mu[i] = viscosity_/d;
output_dmudp[i] = (viscosity_/(d*d))*(1+x) * visc_comp_;
}
} else {
std::fill(output_mu, output_mu + n, viscosity_);
std::fill(output_dmudp, output_dmudp + n, 0.0);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
virtual void B(const int n,
const double* p,
const double* /*z*/,
double* output_B) const
{
if (comp_) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
double x = comp_*(p[i] - ref_press_);
output_B[i] = ref_B_/(1.0 + x + 0.5*x*x);
}
} else {
std::fill(output_B, output_B + n, ref_B_);
}
}
virtual void dBdp(const int n,
const double* p,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
if (comp_) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
double x = comp_*(p[i] - ref_press_);
double d = (1.0 + x + 0.5*x*x);
output_B[i] = ref_B_/d;
output_dBdp[i] = (-ref_B_/(d*d))*(1 + x) * comp_;
}
} else {
std::fill(output_B, output_B + n, ref_B_);
std::fill(output_dBdp, output_dBdp + n, 0.0);
}
}
virtual void b(const int n,
const double* p,
const double* /*r*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
if (comp_) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
double x = comp_*(p[i] - ref_press_);
double d = (1.0 + x + 0.5*x*x);
// b = 1/B = d/ref_B_B;
output_b[i] = d/ref_B_;
output_dbdp[i] = (1 + x) * comp_/ref_B_;
}
} else {
std::fill(output_b, output_b + n, 1/ref_B_);
std::fill(output_dbdp, output_dbdp + n, 0.0);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
virtual void b(const int n,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
if (comp_) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
double x = comp_*(p[i] - ref_press_);
double d = (1.0 + x + 0.5*x*x);
// b = 1/B = d/ref_B_B;
output_b[i] = d/ref_B_;
output_dbdp[i] = (1 + x) * comp_/ref_B_;
}
} else {
std::fill(output_b, output_b + n, 1/ref_B_);
std::fill(output_dbdp, output_dbdp + n, 0.0);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
virtual void rsSat(const int n,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
virtual void rvSat(const int n,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
{
std::fill(output_rvSat, output_rvSat + n, 0.0);
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
virtual void R(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
{
std::fill(output_R, output_R + n, 0.0);
}
virtual void dRdp(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_R,
double* output_dRdp) const
{
std::fill(output_R, output_R + n, 0.0);
std::fill(output_dRdp, output_dRdp + n, 0.0);
}
private:
double ref_press_;
double ref_B_;
double comp_;
double viscosity_;
double visc_comp_;
};
}
#endif // OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED

226
opm/core/props/pvt/SinglePvtDeadSpline.cpp
View File

@ -1,226 +0,0 @@
/*
Copyright 2010, 2011, 2012 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/>.
*/
#include "config.h"
#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/parser/eclipse/Utility/SingleRecordTable.hpp>
#include <algorithm>
// Extra includes for debug dumping of tables.
// #include <boost/lexical_cast.hpp>
// #include <string>
// #include <fstream>
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples)
{
int numRows = pvdoTable.numRows();
// Copy data
const std::vector<double> &press = pvdoTable.getPressureColumn();
const std::vector<double> &B = pvdoTable.getFormationFactorColumn();
const std::vector<double> &visc = pvdoTable.getViscosityColumn();
std::vector<double> B_inv(numRows);
for (int i = 0; i < numRows; ++i) {
B_inv[i] = 1.0 / B[i];
}
buildUniformMonotoneTable(press, B_inv, samples, b_);
buildUniformMonotoneTable(press, visc, samples, viscosity_);
}
SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples)
{
int numRows = pvdgTable.numRows();
// Copy data
const std::vector<double> &press = pvdgTable.getPressureColumn();
const std::vector<double> &B = pvdgTable.getFormationFactorColumn();
const std::vector<double> &visc = pvdgTable.getViscosityColumn();
std::vector<double> B_inv(numRows);
for (int i = 0; i < numRows; ++i) {
B_inv[i] = 1.0 / B[i];
}
buildUniformMonotoneTable(press, B_inv, samples, b_);
buildUniformMonotoneTable(press, visc, samples, viscosity_);
}
// Destructor
SinglePvtDeadSpline::~SinglePvtDeadSpline()
{
}
void SinglePvtDeadSpline::mu(const int n,
const double* p,
const double* /*z*/,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = viscosity_(p[i]);
}
}
void SinglePvtDeadSpline::mu(const int n,
const double* p,
const double* /*r*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = viscosity_(p[i]);
output_dmudp[i] = viscosity_.derivative(p[i]);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void SinglePvtDeadSpline::mu(const int n,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = viscosity_(p[i]);
output_dmudp[i] = viscosity_.derivative(p[i]);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void SinglePvtDeadSpline::B(const int n,
const double* p,
const double* /*z*/,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[i] = 1.0/b_(p[i]);
}
}
void SinglePvtDeadSpline::dBdp(const int n,
const double* p,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
B(n, p, 0, output_B);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
double Bg = output_B[i];
output_dBdp[i] = -Bg*Bg*b_.derivative(p[i]);
}
}
void SinglePvtDeadSpline::b(const int n,
const double* p,
const double* /*r*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_b[i] = b_(p[i]);
output_dbdp[i] = b_.derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void SinglePvtDeadSpline::b(const int n,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_b[i] = b_(p[i]);
output_dbdp[i] = b_.derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void SinglePvtDeadSpline::rsSat(const int n,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
void SinglePvtDeadSpline::rvSat(const int n,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
{
std::fill(output_rvSat, output_rvSat + n, 0.0);
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
void SinglePvtDeadSpline::R(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
{
std::fill(output_R, output_R + n, 0.0);
}
void SinglePvtDeadSpline::dRdp(const int n,
const double* /*p*/,
const double* /*z*/,
double* output_R,
double* output_dRdp) const
{
std::fill(output_R, output_R + n, 0.0);
std::fill(output_dRdp, output_dRdp + n, 0.0);
}
}

48
opm/core/props/pvt/SinglePvtInterface.cpp
View File

@ -1,48 +0,0 @@
/*
Copyright 2010, 2011, 2012 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/>.
*/
#include "config.h"
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
namespace Opm
{
SinglePvtInterface::SinglePvtInterface()
: num_phases_(MaxNumPhases)
{
for (int i = 0; i < MaxNumPhases; ++i) {
phase_pos_[i] = i;
}
}
SinglePvtInterface::~SinglePvtInterface()
{
}
void SinglePvtInterface::setPhaseConfiguration(const int num_phases, const int* phase_pos)
{
num_phases_ = num_phases;
for (int i = 0; i < MaxNumPhases; ++i) {
phase_pos_[i] = phase_pos[i];
}
}
} // namespace Opm

492
opm/core/props/pvt/SinglePvtLiveGas.cpp
View File

@ -1,492 +0,0 @@
//===========================================================================
//
// File: MiscibilityLiveGas.cpp
//
// Created: Wed Feb 10 09:21:53 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
/*
Copyright 2010 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/>.
*/
#include "config.h"
#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <algorithm>
#include <iostream>
namespace Opm
{
using Opm::linearInterpolation;
using Opm::linearInterpolationDerivative;
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
SinglePvtLiveGas::SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable)
{
// GAS, PVTG
saturated_gas_table_.resize(4);
saturated_gas_table_[0] = pvtgTable.getOuterTable()->getPressureColumn();
saturated_gas_table_[1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn();
saturated_gas_table_[2] = pvtgTable.getOuterTable()->getGasViscosityColumn();
saturated_gas_table_[3] = pvtgTable.getOuterTable()->getOilSolubilityColumn();
int sz = pvtgTable.getOuterTable()->numRows();
undersat_gas_tables_.resize(sz);
for (int i=0; i<sz; ++i) {
const auto &undersatTable = *pvtgTable.getInnerTable(i);
undersat_gas_tables_[i].resize(3);
undersat_gas_tables_[i][0] = undersatTable.getOilSolubilityColumn();
undersat_gas_tables_[i][1] = undersatTable.getGasFormationFactorColumn();
undersat_gas_tables_[i][2] = undersatTable.getGasViscosityColumn();
}
// Bg -> 1/Bg
for (int i=0; i<sz; ++i) {
saturated_gas_table_[1][i] = 1.0/saturated_gas_table_[1][i];
for (size_t j=0; j<undersat_gas_tables_[i][1].size(); ++j) {
undersat_gas_tables_[i][1][j] = 1.0/undersat_gas_tables_[i][1][j];
}
}
}
// Destructor
SinglePvtLiveGas::~SinglePvtLiveGas()
{
}
void SinglePvtLiveGas::mu(const int n,
const double* p,
const double* z,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = miscible_gas(p[i], z + num_phases_*i, 2, false);
}
}
/// Viscosity and its derivatives as a function of p and r.
void SinglePvtLiveGas::mu(const int /*n*/,
const double* /*p*/,
const double* /*r*/,
double* /*output_mu*/,
double* /*output_dmudp*/,
double* /*output_dmudr*/) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
/// Viscosity and its derivatives as a function of p and r.
void SinglePvtLiveGas::mu(const int n,
const double* p,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
output_mu[i] = miscible_gas(p[i], r[i], cnd,2, 0);
output_dmudp[i] = miscible_gas(p[i], r[i], cnd, 2, 1);
output_dmudr[i] = miscible_gas(p[i], r[i], cnd, 2, 2);
}
}
/// Formation volume factor as a function of p and z.
void SinglePvtLiveGas::B(const int n,
const double* p,
const double* z,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[i] = evalB(p[i], z + num_phases_*i);
}
}
/// Formation volume factor and p-derivative as functions of p and z.
void SinglePvtLiveGas::dBdp(const int n,
const double* p,
const double* z,
double* output_B,
double* output_dBdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalBDeriv(p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
}
}
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
void SinglePvtLiveGas::b(const int /*n*/,
const double* /*p*/,
const double* /*r*/,
double* /*output_b*/,
double* /*output_dbdp*/,
double* /*output_dbdr*/) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
void SinglePvtLiveGas::b(const int n,
const double* p,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
output_b[i] = miscible_gas(p[i], r[i], cnd, 1, 0);
output_dbdp[i] = miscible_gas(p[i], r[i], cnd, 1, 1);
output_dbdr[i] = miscible_gas(p[i], r[i], cnd, 1, 2);
}
}
/// Gas resolution and its derivatives at bublepoint as a function of p.
void SinglePvtLiveGas::rvSat(const int n,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const
{
for (int i = 0; i < n; ++i) {
output_rvSat[i] = linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[3],p[i]);
output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[0],
saturated_gas_table_[3],p[i]);
}
}
void SinglePvtLiveGas::rsSat(const int n,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
/// Solution factor as a function of p and z.
void SinglePvtLiveGas::R(const int n,
const double* p,
const double* z,
double* output_R) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[i] = evalR(p[i], z + num_phases_*i);
}
}
/// Solution factor and p-derivative as functions of p and z.
void SinglePvtLiveGas::dRdp(const int n,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalRDeriv(p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
}
}
// ---- Private methods ----
double SinglePvtLiveGas::evalB(const double press, const double* surfvol) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
// To handle no-gas case.
return 1.0;
}
return 1.0/miscible_gas(press, surfvol, 1, false);
}
void SinglePvtLiveGas::evalBDeriv(const double press, const double* surfvol,
double& Bval, double& dBdpval) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
// To handle no-gas case.
Bval = 1.0;
dBdpval = 0.0;
return;
}
Bval = evalB(press, surfvol);
dBdpval = -Bval*Bval*miscible_gas(press, surfvol, 1, true);
}
double SinglePvtLiveGas::evalR(const double press, const double* surfvol) const
{
if (surfvol[phase_pos_[Liquid]] == 0.0) {
// To handle no-gas case.
return 0.0;
}
double satR = linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[3], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) {
// Saturated case
return satR;
} else {
// Undersaturated case
return maxR;
}
}
void SinglePvtLiveGas::evalRDeriv(const double press, const double* surfvol,
double& Rval, double& dRdpval) const
{
if (surfvol[phase_pos_[Liquid]] == 0.0) {
// To handle no-gas case.
Rval = 0.0;
dRdpval = 0.0;
return;
}
double satR = linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[3], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) {
// Saturated case
Rval = satR;
dRdpval = linearInterpolationDerivative(saturated_gas_table_[0],
saturated_gas_table_[3],
press);
} else {
// Undersaturated case
Rval = maxR;
dRdpval = 0.0;
}
}
double SinglePvtLiveGas::miscible_gas(const double press,
const double* surfvol,
const int item,
const bool deriv) const
{
int section;
double Rval = linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[3], press,
section);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (deriv) {
if (Rval < maxR ) { // Saturated case
return linearInterpolationDerivative(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
int is = section;
if (undersat_gas_tables_[is][0].size() < 2) {
double val = (saturated_gas_table_[item][is+1]
- saturated_gas_table_[item][is]) /
(saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][is]);
return val;
}
double val1 =
linearInterpolation(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
maxR);
double val2 =
linearInterpolation(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
maxR);
double val = (val2 - val1)/
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
return val;
}
} else {
if (Rval < maxR ) { // Saturated case
return linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
int is = section;
// Extrapolate from first table section
if (is == 0 && press < saturated_gas_table_[0][0]) {
return linearInterpolation(undersat_gas_tables_[0][0],
undersat_gas_tables_[0][item],
maxR);
}
// Extrapolate from last table section
int ltp = saturated_gas_table_[0].size() - 1;
if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
return linearInterpolation(undersat_gas_tables_[ltp][0],
undersat_gas_tables_[ltp][item],
maxR);
}
// Interpolate between table sections
double w = (press - saturated_gas_table_[0][is]) /
(saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][is]);
if (undersat_gas_tables_[is][0].size() < 2) {
double val = saturated_gas_table_[item][is] +
w*(saturated_gas_table_[item][is+1] -
saturated_gas_table_[item][is]);
return val;
}
double val1 =
linearInterpolation(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
maxR);
double val2 =
linearInterpolation(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
maxR);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
double SinglePvtLiveGas::miscible_gas(const double press,
const double r,
const PhasePresence& cond,
const int item,
const int deriv) const
{
const bool isSat = cond.hasFreeOil();
// Derivative w.r.t p
if (deriv == 1) {
if (isSat) { // Saturated case
return linearInterpolationDerivative(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_gas_table_[0], press);
if (undersat_gas_tables_[is][0].size() < 2) {
double val = (saturated_gas_table_[item][is+1]
- saturated_gas_table_[item][is]) /
(saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][is]);
return val;
}
double val1 =
linearInterpolation(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
r);
double val2 =
linearInterpolation(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
r);
double val = (val2 - val1)/
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
return val;
}
} else if (deriv == 2){
if (isSat) {
return 0;
} else {
int is = tableIndex(saturated_gas_table_[0], press);
double w = (press - saturated_gas_table_[0][is]) /
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
assert(undersat_gas_tables_[is][0].size() >= 2);
assert(undersat_gas_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
r);
double val2 =
linearInterpolationDerivative(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
r);
double val = val1 + w * (val2 - val1);
return val;
}
} else {
if (isSat) { // Saturated case
return linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_gas_table_[0], press);
// Extrapolate from first table section
if (is == 0 && press < saturated_gas_table_[0][0]) {
return linearInterpolation(undersat_gas_tables_[0][0],
undersat_gas_tables_[0][item],
r);
}
// Extrapolate from last table section
//int ltp = saturated_gas_table_[0].size() - 1;
//if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
// return linearInterpolation(undersat_gas_tables_[ltp][0],
// undersat_gas_tables_[ltp][item],
// r);
//}
// Interpolate between table sections
double w = (press - saturated_gas_table_[0][is]) /
(saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][is]);
if (undersat_gas_tables_[is][0].size() < 2) {
double val = saturated_gas_table_[item][is] +
w*(saturated_gas_table_[item][is+1] -
saturated_gas_table_[item][is]);
return val;
}
double val1 =
linearInterpolation(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
r);
double val2 =
linearInterpolation(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
r);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

551
opm/core/props/pvt/SinglePvtLiveOil.cpp
View File

@ -1,551 +0,0 @@
/*
Copyright 2010, 2011, 2012 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/>.
*/
#include "config.h"
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <algorithm>
namespace Opm
{
using Opm::linearInterpolation;
using Opm::linearInterpolationDerivative;
using Opm::tableIndex;
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
SinglePvtLiveOil::SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable)
{
const auto saturatedPvto = pvtoTable.getOuterTable();
// OIL, PVTO
saturated_oil_table_.resize(4);
const int sz = saturatedPvto->numRows();
for (int k=0; k<4; ++k) {
saturated_oil_table_[k].resize(sz);
}
for (int i=0; i<sz; ++i) {
saturated_oil_table_[0][i] = saturatedPvto->getPressureColumn()[i]; // p
saturated_oil_table_[1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo
saturated_oil_table_[2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o
saturated_oil_table_[3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs
}
undersat_oil_tables_.resize(sz);
for (int i=0; i<sz; ++i) {
const auto undersaturatedPvto = pvtoTable.getInnerTable(i);
undersat_oil_tables_[i].resize(3);
int tsize = undersaturatedPvto->numRows();
undersat_oil_tables_[i][0].resize(tsize);
undersat_oil_tables_[i][1].resize(tsize);
undersat_oil_tables_[i][2].resize(tsize);
for (int j=0; j<tsize; ++j) {
undersat_oil_tables_[i][0][j] = undersaturatedPvto->getPressureColumn()[j]; // p
undersat_oil_tables_[i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j]; // 1/Bo
undersat_oil_tables_[i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j]; // mu_o
}
}
// Complete undersaturated tables by extrapolating from existing data
// as is done in Eclipse and Mrst
int iNext = -1;
for (int i=0; i<sz; ++i) {
// Skip records already containing undersaturated data
if (undersat_oil_tables_[i][0].size() > 1) {
continue;
}
// Look ahead for next record containing undersaturated data
if (iNext < i) {
iNext = i+1;
while (iNext<sz && undersat_oil_tables_[iNext][0].size() < 2) {
++iNext;
}
if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
}
// Add undersaturated data to current record while maintaining compressibility and viscosibility
typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
for (sz_t j=1; j<undersat_oil_tables_[iNext][0].size(); ++j) {
double diffPressure = undersat_oil_tables_[iNext][0][j]-undersat_oil_tables_[iNext][0][j-1];
double pressure = undersat_oil_tables_[i][0].back()+diffPressure;
undersat_oil_tables_[i][0].push_back(pressure);
double compr = (1.0/undersat_oil_tables_[iNext][1][j]-1.0/undersat_oil_tables_[iNext][1][j-1])
/ (0.5*(1.0/undersat_oil_tables_[iNext][1][j]+1.0/undersat_oil_tables_[iNext][1][j-1]));
double B = (1.0/undersat_oil_tables_[i][1].back())*(1.0+0.5*compr)/(1.0-0.5*compr);
undersat_oil_tables_[i][1].push_back(1.0/B);
double visc = (undersat_oil_tables_[iNext][2][j]-undersat_oil_tables_[iNext][2][j-1])
/ (0.5*(undersat_oil_tables_[iNext][2][j]+undersat_oil_tables_[iNext][2][j-1]));
double mu = (undersat_oil_tables_[i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc);
undersat_oil_tables_[i][2].push_back(mu);
}
}
}
/// Destructor.
SinglePvtLiveOil::~SinglePvtLiveOil()
{
}
/// Viscosity as a function of p and z.
void SinglePvtLiveOil::mu(const int n,
const double* p,
const double* z,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = miscible_oil(p[i], z + num_phases_*i, 2, false);
}
}
/// Viscosity and its derivatives as a function of p and r.
void SinglePvtLiveOil::mu(const int n,
const double* p,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = miscible_oil(p[i], r[i], 2, 0);
output_dmudp[i] = miscible_oil(p[i], r[i], 2, 1);
output_dmudr[i] = miscible_oil(p[i], r[i], 2, 2);
}
}
/// Viscosity and its derivatives as a function of p and r.
void SinglePvtLiveOil::mu(const int n,
const double* p,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
output_mu[i] = miscible_oil(p[i], r[i], cnd, 2, 0);
output_dmudp[i] = miscible_oil(p[i], r[i], cnd, 2, 1);
output_dmudr[i] = miscible_oil(p[i], r[i], cnd, 2, 2);
}
}
/// Formation volume factor as a function of p and z.
void SinglePvtLiveOil::B(const int n,
const double* p,
const double* z,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[i] = evalB(p[i], z + num_phases_*i);
}
}
/// Formation volume factor and p-derivative as functions of p and z.
void SinglePvtLiveOil::dBdp(const int n,
const double* p,
const double* z,
double* output_B,
double* output_dBdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalBDeriv(p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
}
}
void SinglePvtLiveOil::b(const int n,
const double* p,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_b[i] = miscible_oil(p[i], r[i], 1, 0);
output_dbdp[i] = miscible_oil(p[i], r[i], 1, 1);
output_dbdr[i] = miscible_oil(p[i], r[i], 1, 2);
}
}
void SinglePvtLiveOil::b(const int n,
const double* p,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
output_b[i] = miscible_oil(p[i], r[i], cnd, 1, 0);
output_dbdp[i] = miscible_oil(p[i], r[i], cnd, 1, 1);
output_dbdr[i] = miscible_oil(p[i], r[i], cnd, 1, 2);
}
}
void SinglePvtLiveOil::rsSat(const int n,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const
{
for (int i = 0; i < n; ++i) {
output_rsSat[i] = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3],p[i]);
output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[3],p[i]);
}
}
void SinglePvtLiveOil::rvSat(const int n,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
{
std::fill(output_rvSat, output_rvSat + n, 0.0);
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
/// Solution factor as a function of p and z.
void SinglePvtLiveOil::R(const int n,
const double* p,
const double* z,
double* output_R) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[i] = evalR(p[i], z + num_phases_*i);
}
}
/// Solution factor and p-derivative as functions of p and z.
void SinglePvtLiveOil::dRdp(const int n,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalRDeriv(p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
}
}
// ---- Private methods ----
double SinglePvtLiveOil::evalB(double press, const double* surfvol) const
{
// if (surfvol[phase_pos_[Liquid]] == 0.0) return 1.0; // To handle no-oil case.
return 1.0/miscible_oil(press, surfvol, 1, false);
}
void SinglePvtLiveOil::evalBDeriv(const double press, const double* surfvol,
double& Bval, double& dBdpval) const
{
Bval = evalB(press, surfvol);
dBdpval = -Bval*Bval*miscible_oil(press, surfvol, 1, true);
}
double SinglePvtLiveOil::evalR(double press, const double* surfvol) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
return 0.0;
}
double Rval = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) { // Saturated case
return Rval;
} else {
return maxR; // Undersaturated case
}
}
void SinglePvtLiveOil::evalRDeriv(const double press, const double* surfvol,
double& Rval, double& dRdpval) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
Rval = 0.0;
dRdpval = 0.0;
return;
}
Rval = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) {
// Saturated case
dRdpval = linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[3],
press);
} else {
// Undersaturated case
Rval = maxR;
dRdpval = 0.0;
}
}
double SinglePvtLiveOil::miscible_oil(const double press,
const double* surfvol,
const int item,
const bool deriv) const
{
int section;
double Rval = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3],
press, section);
double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (deriv) {
if (Rval < maxR ) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], maxR);
double w = (maxR - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
} else {
if (Rval < maxR ) { // Saturated case
return linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[3], maxR);
double w = (maxR - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
double SinglePvtLiveOil::miscible_oil(const double press,
const double r,
const int item,
const int deriv) const
{
int section;
double Rval = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3],
press, section);
// derivative with respect to frist component (pressure)
if (deriv == 1) {
if (Rval <= r ) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], r);
double w = (r - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
// derivative with respect to second component (r)
} else if (deriv == 2) {
if (Rval <= r ) { // Saturated case
return 0;
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], r);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = (val2 - val1)/(saturated_oil_table_[3][is+1]-saturated_oil_table_[3][is]);
return val;
}
} else {
if (Rval <= r ) { // Saturated case
return linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[3], r);
double w = (r - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
double SinglePvtLiveOil::miscible_oil(const double press,
const double r,
const PhasePresence& cond,
const int item,
const int deriv) const
{
const bool isSat = cond.hasFreeGas();
// derivative with respect to frist component (pressure)
if (deriv == 1) {
if (isSat) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], r);
double w = (r - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
// derivative with respect to second component (r)
} else if (deriv == 2) {
if (isSat) { // Saturated case
return 0;
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], r);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = (val2 - val1)/(saturated_oil_table_[3][is+1]-saturated_oil_table_[3][is]);
return val;
}
} else {
if (isSat) { // Saturated case
return linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[3], r);
double w = (r - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

2
opm/core/simulator/EquilibrationHelpers.hpp
View File

@ -139,7 +139,7 @@ namespace Opm
props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
std::vector<double> rho(np, 0.0);
props_.density(1, &A[0], &rho[0]);
props_.density(1, &A[0], &c_[0], &rho[0]);
return rho;
}

2
opm/core/simulator/initStateEquil.hpp
View File

@ -369,8 +369,6 @@ namespace Opm
const UnstructuredGrid& G ,
const double grav)
{
typedef Miscibility::NoMixing NoMix;
for (typename RMap::RegionId
r = 0, nr = reg.numRegions();
r < nr; ++r)

2
opm/core/simulator/initState_impl.hpp
View File

@ -184,7 +184,7 @@ namespace Opm
double A[4] = { 0.0 };
props_.matrix(1, &pressure, surfvol[phase], cells, A, 0);
double rho[2] = { 0.0 };
props_.density(1, A, rho);
props_.density(1, A, cells, rho);
return rho[phase];
}
};

2
opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp
View File

@ -418,7 +418,7 @@ namespace Opm
assert(np == 2);
const int dim = grid_.dimensions;
density_.resize(nc*np);
props_.density(grid_.number_of_cells, &A_[0], &density_[0]);
props_.density(grid_.number_of_cells, &A_[0], /*cellIndices=*/NULL, &density_[0]);
std::fill(gravflux_.begin(), gravflux_.end(), 0.0);
for (int f = 0; f < nf; ++f) {
const int* c = &grid_.face_cells[2*f];

3
opm/core/utility/miscUtilitiesBlackoil.cpp
View File

@ -285,7 +285,8 @@ namespace Opm
/// @brief Computes saturation from surface volume
void computeSaturation(const BlackoilPropertiesInterface& props,
BlackoilState& state){
BlackoilState& state)
{
const int np = props.numPhases();
const int nc = props.numCells();

4
tests/capillary.DATA
View File

@ -19,6 +19,10 @@ PVDG
200 0.005 0.2
/
PVTW
1.0 1.0 4.0E-5 0.96 0.0
/
SWOF
0.2 0 1 0.4
1 1 0 0.1

4
tests/deadfluids.DATA
View File

@ -19,6 +19,10 @@ PVDG
200 0.02 0.2
/
PVTW
1.0 1.0 4.0E-5 0.96 0.0
/
SWOF
0 0 1 0
1 1 0 0

13
tests/liveoil.DATA
View File

@ -9,6 +9,19 @@ WATER
FIELD
-- tests for the PVT functions need a grid because the OPM-API for the
-- PVT functions assumes the presence of compressed grid cells,
-- i.e. we need to be able to map from compressed to logical cartesian
-- cell indices to find out the PVT region of a cell
DIMENS
3 3 3 /
DXV
1 2 3 /
DYV
1 2 3 /
DZV
1 2 3 /
PVTO
-- Rs Pbub Bo Vo
.0 14.7 1.0000 1.20 /

111
tests/test_blackoilfluid.cpp
View File

@ -1,15 +1,17 @@
#include <config.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/props/pvt/SinglePvtConstCompr.hpp>
#include <opm/core/props/pvt/SinglePvtDead.hpp>
#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
#include <opm/core/props/pvt/PvtConstCompr.hpp>
#include <opm/core/props/pvt/PvtDead.hpp>
#include <opm/core/props/pvt/PvtDeadSpline.hpp>
#include <opm/core/props/pvt/PvtLiveOil.hpp>
#include <opm/core/props/pvt/PvtLiveGas.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/grid.h>
#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
#include <opm/parser/eclipse/Utility/PvtgTable.hpp>
@ -36,63 +38,68 @@ using namespace Opm;
using namespace std;
std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr deck,PhaseUsage phase_usage_){
std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, PhaseUsage phase_usage_){
Opm::GridManager grid(deck);
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
int samples = 0;
std::vector<std::shared_ptr<SinglePvtInterface> > props_;
std::vector<std::shared_ptr<PvtInterface> > props_;
// Set the properties.
props_.resize(phase_usage_.num_phases);
// Water PVT
if (phase_usage_.phase_used[Aqua]) {
if (deck->hasKeyword("PVTW")) {
Opm::PvtwTable pvtwTable(deck->getKeyword("PVTW"));
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable));
std::shared_ptr<PvtConstCompr> pvtw(new PvtConstCompr);
pvtw->initFromWater(deck->getKeyword("PVTW"));
props_[phase_usage_.phase_pos[Aqua]] = pvtw;
} else {
// Eclipse 100 default.
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
props_[phase_usage_.phase_pos[Aqua]].reset(new PvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
}
}
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
if (deck->hasKeyword("PVDO")) {
Opm::PvdoTable pvdoTable(deck->getKeyword("PVDO"));
Opm::DeckKeywordConstPtr pvdoKeyword(deck->getKeyword("PVDO"));
if (samples > 0) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromOil(pvdoKeyword, samples);
props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
} else {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable));
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromOil(pvdoKeyword);
props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
}
} else if (deck->hasKeyword("PVTO")) {
Opm::PvtoTable pvtoTable(deck->getKeyword("PVTO"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable));
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO")));
} else if (deck->hasKeyword("PVCDO")) {
Opm::PvcdoTable pvcdoTable(deck->getKeyword("PVCDO"), /*tableIdx=*/0);
std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
pvcdo->initFromOil(deck->getKeyword("PVCDO"));
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable));
props_[phase_usage_.phase_pos[Liquid]] = pvcdo;
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n");
}
}
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
if (deck->hasKeyword("PVDG")) {
Opm::PvdgTable pvdgTable(deck->getKeyword("PVDG"));
Opm::DeckKeywordConstPtr pvdgKeyword(deck->getKeyword("PVDG"));
if (samples > 0) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromGas(pvdgKeyword, samples);
props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
} else {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable));
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromGas(pvdgKeyword);
props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
}
} else if (deck->hasKeyword("PVTG")) {
Opm::PvtgTable pvtgTable(deck->getKeyword("PVTG"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable));
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG")));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
}
@ -101,8 +108,8 @@ std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr dec
return props_;
}
void testmu(const double reltol, int n, int np, std::vector<double> p, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<SinglePvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
void testmu(const double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{
std::vector<double> mu(n);
std::vector<double> dmudp(n);
@ -115,8 +122,8 @@ void testmu(const double reltol, int n, int np, std::vector<double> p, std::vect
// test mu
for (int phase = 0; phase < np; ++phase) {
props_[phase]->mu(n, &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
props_[phase]->mu(n, &p[0], &z[0], &mu[0]);
props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &z[0], &mu[0]);
dmudp_diff = (mu_new[1]-mu_new[0])/(p[1]-p[0]);
dmudr_diff = (mu_new[2]-mu_new[0])/(r[2]-r[0]);
dmudp_diff_u = (mu_new[4]-mu_new[3])/(p[4]-p[3]);
@ -137,8 +144,8 @@ void testmu(const double reltol, int n, int np, std::vector<double> p, std::vect
}
}
void testb(const double reltol, int n, int np, std::vector<double> p, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<SinglePvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
void testb(const double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{
// test b
std::vector<double> b(n);
@ -154,8 +161,8 @@ void testb(const double reltol, int n, int np, std::vector<double> p, std::vecto
double dbdr_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase]->b(n, &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
props_[phase]->dBdp(n, &p[0], &z[0], &B[0], &dBdp[0]);
props_[phase]->b(n, &pvtTableIdx[0], &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
props_[phase]->dBdp(n, &pvtTableIdx[0], &p[0], &z[0], &B[0], &dBdp[0]);
dbdp_diff = (b[1]-b[0])/(p[1]-p[0]);
dbdr_diff = (b[2]-b[0])/(r[2]-r[0]);
dbdp_diff_u = (b[4]-b[3])/(p[4]-p[3]);
@ -180,7 +187,7 @@ void testb(const double reltol, int n, int np, std::vector<double> p, std::vecto
}
}
void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
void testrsSat(double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<std::shared_ptr<PvtInterface> > props_){
// test bublepoint pressure
std::vector<double> rs(n);
std::vector<double> drsdp(n);
@ -188,7 +195,7 @@ void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<
double drsdp_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase] ->rsSat(n, &p[0], &rs[0], &drsdp[0]);
props_[phase] ->rsSat(n, &pvtTableIdx[0], &p[0], &rs[0], &drsdp[0]);
drsdp_diff = (rs[1]-rs[0])/(p[1]-p[0]);
drsdp_diff_u = (rs[4]-rs[3])/(p[4]-p[3]);
@ -202,7 +209,7 @@ void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<
}
}
void testrvSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
void testrvSat(double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<std::shared_ptr<PvtInterface> > props_){
// test rv saturated
std::vector<double> rv(n);
std::vector<double> drvdp(n);
@ -210,7 +217,7 @@ void testrvSat(double reltol, int n, int np, std::vector<double> p, std::vector<
double drvdp_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase] ->rvSat(n, &p[0], &rv[0], &drvdp[0]);
props_[phase] ->rvSat(n, &pvtTableIdx[0], &p[0], &rv[0], &drvdp[0]);
drvdp_diff = (rv[1]-rv[0])/(p[1]-p[0]);
drvdp_diff_u = (rv[4]-rv[3])/(p[4]-p[3]);
@ -234,15 +241,18 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr deck(parser->parseFile(filename));
// setup pvt interface
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(deck,phase_usage_);
std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck,phase_usage_);
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
// approximation of the derivatives.
const int n = 6;
const int np = phase_usage_.num_phases;
std::vector<int> pvtRegionIdx(n, 0);
// the tolerance for acceptable difference in values
const double reltol = 1e-9;
@ -287,13 +297,13 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
}
testmu(reltol, n, np, p, r,z, props_, condition);
testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition);
testb(reltol,n,np,p,r,z,props_,condition);
testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition);
testrsSat(reltol,n,np,p,props_);
testrsSat(reltol,n,np,pvtRegionIdx,p,props_);
testrvSat(reltol,n,np,p,props_);
testrvSat(reltol,n,np,pvtRegionIdx,p,props_);
}
@ -311,13 +321,14 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
// setup pvt interface
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(deck,phase_usage_);
std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck,phase_usage_);
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
// approximation of the derivatives.
const int n = 6;
const int np = phase_usage_.num_phases;
std::vector<int> pvtRegionIdx(n, 0);
// the tolerance for acceptable difference in values
const double reltol = 1e-9;
@ -362,12 +373,12 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
}
testmu(reltol, n, np, p, r,z, props_, condition);
testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition);
testb(reltol,n,np,p,r,z,props_,condition);
testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition);
testrsSat(reltol,n,np,p,props_);
testrsSat(reltol,n,np,pvtRegionIdx,p,props_);
testrvSat(reltol,n,np,p,props_);
testrvSat(reltol,n,np,pvtRegionIdx,p,props_);
}

8
tests/test_equil.cpp
View File

@ -428,13 +428,13 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
// the true answer or something else.
const double reltol = 1.0e-6;
BOOST_CHECK_CLOSE(pressures[0][first] , 1.469769063e7 , reltol);
BOOST_CHECK_CLOSE(pressures[0][last ] , 1.545e7 , reltol);
BOOST_CHECK_CLOSE(pressures[1][last] , 1.546e7 , reltol);
BOOST_CHECK_CLOSE(pressures[0][last ] , 15452880.328284413 , reltol);
BOOST_CHECK_CLOSE(pressures[1][last] , 15462880.328284413 , reltol);
const auto& sats = comp.saturation();
const std::vector<double> s[3]{
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.425893333333, 0.774026666666, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0.00736, 0.792746666666, 0.8, 0.8, 0.8, 0.8, 0.574106666666, 0.225973333333, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.42192000000000002, 0.77802666666666664, 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.57807999999999993, 0.22197333333333336, 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) {

13
tests/wetgas.DATA
View File

@ -9,6 +9,19 @@ GAS
FIELD
-- tests for the PVT functions need a grid because the OPM-API for the
-- PVT functions assumes the presence of compressed grid cells,
-- i.e. we need to be able to map from compressed to logical cartesian
-- cell indices to find out the PVT region of a cell
DIMENS
3 3 3 /
DXV
1 2 3 /
DYV
1 2 3 /
DZV
1 2 3 /
-- PVT PROPERTIES OF DRY GAS (NO VAPOURISED OIL)
-- FROM SPE3 Blackoil Kleppe
--