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

Merge pull request #657 from GitPaean/PVT_viscosity_try

Browse Source 
Obtaining viscosity by doing a linear interpolation of the inverse of B*mu instead of mu directly.
This commit is contained in:
Atgeirr Flø Rasmussen 2014-09-23 22:44:37 +02:00
commit c3551dfdad
5 changed files with 215 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
opm/core/props/pvt/PvtDead.cpp
View File

@ -41,6 +41,7 @@ namespace Opm
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
inverseBmu_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::PvdoTable& pvdoTable = pvdoTables[regionIdx];
@ -51,12 +52,19 @@ namespace Opm
const std::vector<double>& visc = pvdoTable.getViscosityColumn();
const int sz = b.size();
std::vector<double> bInv(sz);
std::vector<double> inverseB(sz);
for (int i = 0; i < sz; ++i) {
bInv[i] = 1.0 / b[i];
inverseB[i] = 1.0 / b[i];
}
b_[regionIdx] = NonuniformTableLinear<double>(press, bInv);
std::vector<double> inverseBmu(sz);
for (int i = 0; i < sz; ++i) {
inverseBmu[i] = 1.0 / (b[i] * visc[i]);
}
b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu);
}
}
@ -68,6 +76,7 @@ namespace Opm
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
inverseBmu_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::PvdgTable& pvdgTable = pvdgTables[regionIdx];
@ -78,12 +87,19 @@ namespace Opm
const std::vector<double>& visc = pvdgTable.getViscosityColumn();
const int sz = b.size();
std::vector<double> bInv(sz);
std::vector<double> inverseB(sz);
for (int i = 0; i < sz; ++i) {
bInv[i] = 1.0 / b[i];
inverseB[i] = 1.0 / b[i];
}
b_[regionIdx] = NonuniformTableLinear<double>(press, bInv);
std::vector<double> inverseBmu(sz);
for (int i = 0; i < sz; ++i) {
inverseBmu[i] = 1.0 / (b[i] * visc[i]);
}
b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu);
}
}
@ -103,7 +119,9 @@ namespace Opm
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
double tempInvB = b_[regionIdx](p[i]);
double tempInvBmu = inverseBmu_[regionIdx](p[i]);
output_mu[i] = tempInvB / tempInvBmu;
}
}
@ -118,8 +136,11 @@ namespace Opm
// #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]);
double tempInvB = b_[regionIdx](p[i]);
double tempInvBmu = inverseBmu_[regionIdx](p[i]);
output_mu[i] = tempInvB / tempInvBmu;
output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i])
- tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
@ -137,8 +158,12 @@ namespace Opm
// #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]);
double tempInvB = b_[regionIdx](p[i]);
double tempInvBmu = inverseBmu_[regionIdx](p[i]);
output_mu[i] = tempInvB / tempInvBmu;
output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i])
- tempInvB * inverseBmu_[regionIdx].derivative(p[i]))
/ (tempInvBmu * tempInvBmu);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);

1
opm/core/props/pvt/PvtDead.hpp
View File

@ -155,6 +155,7 @@ namespace Opm
// table per PVT region.
std::vector<NonuniformTableLinear<double> > b_;
std::vector<NonuniformTableLinear<double> > viscosity_;
std::vector<NonuniformTableLinear<double> > inverseBmu_;
};
}

72
opm/core/props/pvt/PvtLiveGas.cpp
View File

@ -56,28 +56,40 @@ namespace Opm
const Opm::PvtgTable& pvtgTable = pvtgTables[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();
// saturated_gas_table_[pvtTableIdx].resize(4);
// Adding one extra line to PVTG to store 1./(Bg*mu_g)
saturated_gas_table_[pvtTableIdx].resize(5);
saturated_gas_table_[pvtTableIdx][0] = pvtgTable.getOuterTable()->getPressureColumn(); // Pressure
saturated_gas_table_[pvtTableIdx][1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); // Bg
saturated_gas_table_[pvtTableIdx][2] = pvtgTable.getOuterTable()->getGasViscosityColumn(); // mu_g
// The number of the columns
int nRows = saturated_gas_table_[pvtTableIdx][2].size();
saturated_gas_table_[pvtTableIdx][3].resize(nRows); // allocate memory for 1/(Bg*mu_g)
saturated_gas_table_[pvtTableIdx][4] = pvtgTable.getOuterTable()->getOilSolubilityColumn(); // Rv
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();
// undersat_gas_tables_[pvtTableIdx][i].resize(3);
undersat_gas_tables_[pvtTableIdx][i].resize(4);
undersat_gas_tables_[pvtTableIdx][i][0] = undersatTable.getOilSolubilityColumn(); // Rv
undersat_gas_tables_[pvtTableIdx][i][1] = undersatTable.getGasFormationFactorColumn(); // Bg
undersat_gas_tables_[pvtTableIdx][i][2] = undersatTable.getGasViscosityColumn(); // mu_g
int nRows = undersat_gas_tables_[pvtTableIdx][i][2].size();
undersat_gas_tables_[pvtTableIdx][i][3].resize(nRows); // allocate memory for 1/(Bg*mu_g)
}
// Bg -> 1/Bg
for (int i=0; i<sz; ++i) {
saturated_gas_table_[pvtTableIdx][1][i] = 1.0/saturated_gas_table_[pvtTableIdx][1][i];
saturated_gas_table_[pvtTableIdx][3][i] = 1.0 / (saturated_gas_table_[pvtTableIdx][1][i]
* saturated_gas_table_[pvtTableIdx][2][i]); // 1/(Bg*mu_g)
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];
undersat_gas_tables_[pvtTableIdx][i][3][j] = 1.0 / (undersat_gas_tables_[pvtTableIdx][i][1][j]
* undersat_gas_tables_[pvtTableIdx][i][2][j]); // 1/(Bg*mu_g)
undersat_gas_tables_[pvtTableIdx][i][1][j] = 1.0 / undersat_gas_tables_[pvtTableIdx][i][1][j];
}
}
}
@ -97,7 +109,10 @@ namespace Opm
{
// #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);
double inverseB = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 1, false);
double inverseBMu = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 3, false);
output_mu[i] = inverseB / inverseBMu;
}
}
@ -126,9 +141,24 @@ namespace Opm
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);
double inverseBMu = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 0);
double inverseB = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 0);
output_mu[i] = inverseB / inverseBMu;
double dinverseBmudp = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 1);
double dinverseBdp = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 1);
output_dmudp[i] = (inverseBMu * dinverseBdp - inverseB * dinverseBmudp)
/ (inverseBMu * inverseBMu);
double dinverseBmudr = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 2);
double dinverseBdr = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 2);
output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
/ (inverseBMu * inverseBMu);
}
}
@ -209,9 +239,9 @@ namespace Opm
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]);
saturated_gas_table_[pvtTableIdx][4],p[i]);
output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][3],p[i]);
saturated_gas_table_[pvtTableIdx][4],p[i]);
}
}
@ -287,7 +317,7 @@ namespace Opm
return 0.0;
}
double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][3], press);
saturated_gas_table_[pvtTableIdx][4], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) {
// Saturated case
@ -308,13 +338,13 @@ namespace Opm
return;
}
double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][3], press);
saturated_gas_table_[pvtTableIdx][4], 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],
saturated_gas_table_[pvtTableIdx][4],
press);
} else {
// Undersaturated case
@ -336,7 +366,7 @@ namespace Opm
int section;
double Rval = linearInterpolation(saturatedGasTable[0],
saturatedGasTable[3], press,
saturatedGasTable[4], press,
section);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (deriv) {

171
opm/core/props/pvt/PvtLiveOil.cpp
View File

@ -46,36 +46,47 @@ namespace Opm
const auto saturatedPvto = pvtoTable.getOuterTable();
// OIL, PVTO
saturated_oil_table_[pvtTableIdx].resize(4);
// saturated_oil_table_[pvtTableIdx].resize(4);
// adding a extra colummn to the PVTO to store 1/(B*mu)
saturated_oil_table_[pvtTableIdx].resize(5);
const int sz = saturatedPvto->numRows();
for (int k=0; k<4; ++k) {
// for (int k=0; k<4; ++k) {
for (int k=0; k<5; ++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
saturated_oil_table_[pvtTableIdx][3][i] = 1.0 / (saturatedPvto->getOilFormationFactorColumn()[i]
* saturatedPvto->getOilViscosityColumn()[i]); // 1/(Bo*mu_o)
saturated_oil_table_[pvtTableIdx][4][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);
// undersat_oil_tables_[pvtTableIdx][i].resize(3);
// adding a extra colummn to the PVTO to store 1/(B*mu)
undersat_oil_tables_[pvtTableIdx][i].resize(4);
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);
undersat_oil_tables_[pvtTableIdx][i][3].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
undersat_oil_tables_[pvtTableIdx][i][3][j] = 1.0 / (undersaturatedPvto->getOilFormationFactorColumn()[j] *
undersaturatedPvto->getOilViscosityColumn()[j]); // 1/(Bo*mu_o)
}
}
// Complete undersaturated tables by extrapolating from existing data
// as is done in Eclipse and Mrst
// TODO: check if the following formulations applying to 1/(Bo*mu_o)
int iNext = -1;
for (int i=0; i<sz; ++i) {
// Skip records already containing undersaturated data
@ -91,6 +102,7 @@ namespace Opm
if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
}
// Add undersaturated data to current record while maintaining compressibility and viscosibility
// TODO: How to add 1/(B*mu) in this way?
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];
@ -104,6 +116,11 @@ namespace Opm
/ (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);
// A try to expolate the 1/BMu with the expolated mu and B
double inverseBMu = 1.0 / (B*mu);
undersat_oil_tables_[pvtTableIdx][i][3].push_back(inverseBMu);
}
}
}
@ -126,7 +143,10 @@ namespace Opm
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);
double inverseB = miscible_oil(p[i], z + num_phases_*i, tableIdx, 1, false);
double inverseBMu = miscible_oil(p[i], z + num_phases_*i, tableIdx, 3, false);
output_mu[i] = inverseB / inverseBMu;
}
}
@ -143,9 +163,23 @@ namespace Opm
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);
double inverseBMu = miscible_oil(p[i], r[i], tableIdx, 3, 0);
double inverseB = miscible_oil(p[i], r[i], tableIdx, 1, 0);
output_mu[i] = inverseB / inverseBMu;
double dinverseBmudp = miscible_oil(p[i], r[i], tableIdx, 3, 1);
double dinverseBdp = miscible_oil(p[i], r[i], tableIdx, 1, 1);
output_dmudp[i] = (inverseBMu * dinverseBdp - inverseB * dinverseBmudp)
/ (inverseBMu * inverseBMu);
double dinverseBmudr = miscible_oil(p[i], r[i], tableIdx, 3, 2);
double dinverseBdr = miscible_oil(p[i], r[i], tableIdx, 1, 2);
output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
/ (inverseBMu * inverseBMu);
}
}
@ -165,9 +199,23 @@ namespace Opm
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);
double inverseBMu = miscible_oil(p[i], r[i], cnd, tableIdx, 3, 0);
double inverseB = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 0);
output_mu[i] = inverseB / inverseBMu;
double dinverseBmudp = miscible_oil(p[i], r[i], cnd, tableIdx, 3, 1);
double dinverseBdp = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 1);
output_dmudp[i] = (inverseBMu * dinverseBdp - inverseB * dinverseBmudp)
/ (inverseBMu * inverseBMu);
double dinverseBmudr = miscible_oil(p[i], r[i], cnd, tableIdx, 3, 2);
double dinverseBdr = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 2);
output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
/ (inverseBMu * inverseBMu);
}
}
@ -258,9 +306,9 @@ namespace Opm
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]);
saturated_oil_table_[tableIdx][4],p[i]);
output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[tableIdx][0],
saturated_oil_table_[tableIdx][3],p[i]);
saturated_oil_table_[tableIdx][4],p[i]);
}
}
@ -334,7 +382,7 @@ namespace Opm
return 0.0;
}
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][3], press);
saturated_oil_table_[pvtTableIdx][4], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) { // Saturated case
return Rval;
@ -353,12 +401,12 @@ namespace Opm
return;
}
Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][3], press);
saturated_oil_table_[pvtTableIdx][4], 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],
saturated_oil_table_[pvtTableIdx][4],
press);
} else {
// Undersaturated case
@ -368,6 +416,7 @@ namespace Opm
}
// TODO: Check if this function need to be adapted to the 1/(B*mu) interpolation.
double PvtLiveOil::miscible_oil(const double press,
const double* surfvol,
const int pvtTableIdx,
@ -376,7 +425,7 @@ namespace Opm
{
int section;
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][3],
saturated_oil_table_[pvtTableIdx][4],
press, section);
double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (deriv) {
@ -385,9 +434,9 @@ namespace Opm
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]);
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], maxR);
double w = (maxR - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
@ -408,9 +457,9 @@ namespace Opm
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]);
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], maxR);
double w = (maxR - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
@ -435,7 +484,7 @@ namespace Opm
{
int section;
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][3],
saturated_oil_table_[pvtTableIdx][4],
press, section);
// derivative with respect to frist component (pressure)
if (deriv == 1) {
@ -444,9 +493,9 @@ namespace Opm
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]);
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
@ -465,7 +514,7 @@ namespace Opm
if (Rval <= r ) { // Saturated case
return 0;
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
@ -477,7 +526,7 @@ namespace Opm
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]);
double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][4][is+1]-saturated_oil_table_[pvtTableIdx][4][is]);
return val;
}
@ -489,9 +538,9 @@ namespace Opm
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]);
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
@ -524,9 +573,9 @@ namespace Opm
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]);
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
@ -545,7 +594,7 @@ namespace Opm
if (isSat) { // Saturated case
return 0;
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
@ -557,35 +606,37 @@ namespace Opm
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]);
double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][4][is+1]-saturated_oil_table_[pvtTableIdx][4][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;
}
// no derivative
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][4], r);
double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][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

28
tests/test_blackoilfluid.cpp
View File

@ -247,8 +247,10 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
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;
// the relative tolerance in percentage for acceptable difference in values
const double reltolper = 1e-9;
// the relative tolerance in percentage for acceptable difference in values for viscosity
const double reltolpermu = 1e-1;
std::vector<double> p(n);
std::vector<double> r(n);
@ -290,13 +292,13 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
}
testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition);
testmu(reltolpermu, n, np, pvtRegionIdx, p, r,z, props_, condition);
testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition);
testb(reltolper,n,np,pvtRegionIdx,p,r,z,props_,condition);
testrsSat(reltol,n,np,pvtRegionIdx,p,props_);
testrsSat(reltolper,n,np,pvtRegionIdx,p,props_);
testrvSat(reltol,n,np,pvtRegionIdx,p,props_);
testrvSat(reltolper,n,np,pvtRegionIdx,p,props_);
}
@ -324,8 +326,10 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
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;
// the relative tolerance in percentage for acceptable difference in values
const double reltolper = 1e-9;
// the relative tolerance in percentage for acceptable difference in values for viscosity
const double reltolpermu = 1e-1;
std::vector<double> p(n);
std::vector<double> r(n);
@ -367,12 +371,12 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
}
testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition);
testmu(reltolpermu, n, np, pvtRegionIdx, p, r,z, props_, condition);
testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition);
testb(reltolper,n,np,pvtRegionIdx,p,r,z,props_,condition);
testrsSat(reltol,n,np,pvtRegionIdx,p,props_);
testrsSat(reltolper,n,np,pvtRegionIdx,p,props_);
testrvSat(reltol,n,np,pvtRegionIdx,p,props_);
testrvSat(reltolper,n,np,pvtRegionIdx,p,props_);
}