598c7a0cb8
The pvt interface is extended to handle wet-gas systems: 1. rvSat is added as a function in the PVT interface 2. SinglePvtLiveGas computes the pvt values and its derivatives 3. The old rbub variable is changed to rsSat for clearity 4. The new interface is tested in test_blackoilfluid with data from liveoil.DATA and wetgas.DATA
355 lines
10 KiB
C++
355 lines
10 KiB
C++
#include <config.h>
|
|
|
|
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
|
|
#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/phaseUsageFromDeck.hpp>
|
|
#include <opm/core/props/BlackoilPhases.hpp>
|
|
#include <opm/core/utility/Units.hpp>
|
|
#include <opm/core/utility/ErrorMacros.hpp>
|
|
|
|
#if HAVE_DYNAMIC_BOOST_TEST
|
|
#define BOOST_TEST_DYN_LINK
|
|
#endif
|
|
#define NVERBOSE // to suppress our messages when throwing
|
|
#define BOOST_TEST_MODULE BlackoilFluidTest
|
|
#define BOOST_TEST_MAIN
|
|
#include <boost/test/unit_test.hpp>
|
|
#include <memory>
|
|
#include <iostream>
|
|
#include <iterator>
|
|
#include <vector>
|
|
#include <string>
|
|
|
|
using namespace Opm;
|
|
using namespace std;
|
|
|
|
|
|
std::vector<std::shared_ptr<SinglePvtInterface> > getProps(const EclipseGridParser deck,PhaseUsage phase_usage_){
|
|
|
|
|
|
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
|
|
int samples = 0;
|
|
|
|
|
|
std::vector<std::shared_ptr<SinglePvtInterface> > props_;
|
|
// Set the properties.
|
|
props_.resize(phase_usage_.num_phases);
|
|
|
|
// Water PVT
|
|
if (phase_usage_.phase_used[Aqua]) {
|
|
if (deck.hasField("PVTW")) {
|
|
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(deck.getPVTW().pvtw_));
|
|
} else {
|
|
// Eclipse 100 default.
|
|
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
|
|
}
|
|
}
|
|
|
|
// Oil PVT
|
|
if (phase_usage_.phase_used[Liquid]) {
|
|
if (deck.hasField("PVDO")) {
|
|
if (samples > 0) {
|
|
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(deck.getPVDO().pvdo_, samples));
|
|
} else {
|
|
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(deck.getPVDO().pvdo_));
|
|
}
|
|
} else if (deck.hasField("PVTO")) {
|
|
|
|
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(deck.getPVTO().pvto_));
|
|
} else if (deck.hasField("PVCDO")) {
|
|
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(deck.getPVCDO().pvcdo_));
|
|
} else {
|
|
OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
|
|
}
|
|
}
|
|
// Gas PVT
|
|
if (phase_usage_.phase_used[Vapour]) {
|
|
if (deck.hasField("PVDG")) {
|
|
if (samples > 0) {
|
|
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(deck.getPVDG().pvdg_, samples));
|
|
} else {
|
|
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(deck.getPVDG().pvdg_));
|
|
}
|
|
} else if (deck.hasField("PVTG")) {
|
|
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(deck.getPVTG().pvtg_));
|
|
} else {
|
|
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
|
|
}
|
|
}
|
|
|
|
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)
|
|
{
|
|
std::vector<double> mu(n);
|
|
std::vector<double> dmudp(n);
|
|
std::vector<double> dmudr(n);
|
|
std::vector<double> mu_new(n);
|
|
double dmudp_diff;
|
|
double dmudr_diff;
|
|
double dmudp_diff_u;
|
|
double dmudr_diff_u;
|
|
|
|
// 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]);
|
|
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]);
|
|
dmudr_diff_u = (mu_new[5]-mu_new[3])/(r[5]-r[3]);
|
|
|
|
for (int i = 0; i < n; ++i){
|
|
BOOST_CHECK_CLOSE(mu_new[i],mu[i],reltol);
|
|
}
|
|
|
|
// saturated case
|
|
BOOST_CHECK_CLOSE(dmudp_diff,dmudp[0],reltol);
|
|
BOOST_CHECK_CLOSE(dmudr_diff,dmudr[0],reltol);
|
|
|
|
// unsaturated case
|
|
BOOST_CHECK_CLOSE(dmudp_diff_u,dmudp[3],reltol);
|
|
BOOST_CHECK_CLOSE(dmudr_diff_u,dmudr[3],reltol);
|
|
|
|
}
|
|
}
|
|
|
|
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)
|
|
{
|
|
// test b
|
|
std::vector<double> b(n);
|
|
std::vector<double> B(n);
|
|
std::vector<double> invB(n);
|
|
std::vector<double> dinvBdp(n);
|
|
std::vector<double> dBdp(n);
|
|
std::vector<double> dbdr(n);
|
|
std::vector<double> dbdp(n);
|
|
double dbdp_diff;
|
|
double dbdr_diff;
|
|
double dbdp_diff_u;
|
|
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]);
|
|
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]);
|
|
dbdr_diff_u = (b[5]-b[3])/(r[5]-r[3]);
|
|
for (int i = 0; i < n; ++i){
|
|
invB[i] = 1/B[i];
|
|
dinvBdp[i] = -1/pow(B[i],2) * dBdp[i];
|
|
}
|
|
|
|
for (int i = 0; i < n; ++i){
|
|
BOOST_CHECK_CLOSE(invB[i],b[i] , reltol);
|
|
BOOST_CHECK_CLOSE(dinvBdp[i],dbdp[i] , reltol);
|
|
|
|
}
|
|
// saturated case
|
|
BOOST_CHECK_CLOSE(dbdp_diff,dbdp[0], reltol);
|
|
BOOST_CHECK_CLOSE(dbdr_diff,dbdr[0], reltol);
|
|
|
|
// unsaturated case
|
|
BOOST_CHECK_CLOSE(dbdp_diff_u,dbdp[3], reltol);
|
|
BOOST_CHECK_CLOSE(dbdr_diff_u,dbdr[3], reltol);
|
|
}
|
|
}
|
|
|
|
void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
|
|
// test bublepoint pressure
|
|
std::vector<double> rs(n);
|
|
std::vector<double> drsdp(n);
|
|
double drsdp_diff;
|
|
double drsdp_diff_u;
|
|
|
|
for (int phase = 0; phase < np; ++phase) {
|
|
props_[phase] ->rsSat(n, &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]);
|
|
|
|
// saturated case
|
|
BOOST_CHECK_CLOSE(drsdp_diff,drsdp[0], reltol);
|
|
|
|
// unsaturad case
|
|
BOOST_CHECK_CLOSE(drsdp_diff_u,drsdp[3], reltol);
|
|
|
|
}
|
|
}
|
|
|
|
void testrvSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
|
|
// test rv saturated
|
|
std::vector<double> rv(n);
|
|
std::vector<double> drvdp(n);
|
|
double drvdp_diff;
|
|
double drvdp_diff_u;
|
|
|
|
for (int phase = 0; phase < np; ++phase) {
|
|
props_[phase] ->rvSat(n, &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]);
|
|
|
|
// saturated case
|
|
BOOST_CHECK_CLOSE(drvdp_diff,drvdp[0], reltol);
|
|
|
|
// unsaturad case
|
|
BOOST_CHECK_CLOSE(drvdp_diff_u,drvdp[3], reltol);
|
|
|
|
}
|
|
}
|
|
|
|
BOOST_AUTO_TEST_CASE(test_liveoil)
|
|
{
|
|
|
|
|
|
// read eclipse deck
|
|
const string filename = "liveoil.DATA";
|
|
cout << "Reading deck: " << filename << endl;
|
|
const EclipseGridParser deck (filename);
|
|
|
|
// setup pvt interface
|
|
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
|
|
std::vector<std::shared_ptr<SinglePvtInterface> > 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;
|
|
|
|
// the tolerance for acceptable difference in values
|
|
const double reltol = 1e-9;
|
|
|
|
std::vector<double> p(n);
|
|
std::vector<double> r(n);
|
|
std::vector<PhasePresence> condition(n);
|
|
std::vector<double> z(n * np);
|
|
|
|
|
|
// Used for forward difference calculations
|
|
const double h_p = 1e4;
|
|
const double h_r = 1;
|
|
|
|
// saturated
|
|
p[0] = 1e7;
|
|
p[1] = p[0] + h_p;
|
|
p[2] = p[0];
|
|
|
|
r[0] = 200;
|
|
r[1] = 200;
|
|
r[2] = 200 + h_r;
|
|
|
|
condition[0].setFreeGas();
|
|
condition[1].setFreeGas();
|
|
condition[2].setFreeGas();
|
|
|
|
|
|
// undersaturated
|
|
p[3] = p[0];
|
|
p[4] = p[1];
|
|
p[5] = p[2];
|
|
|
|
r[3] = 50;
|
|
r[4] = 50;
|
|
r[5] = 50 +h_r;
|
|
|
|
// Corresponing z factors, used to compare with the [p,z] interface
|
|
for (int i = 0; i < n; ++i) {
|
|
z[0+i*np] = 0; z[1+i*np] = 1;
|
|
z[2+i*np] = r[i];
|
|
|
|
}
|
|
|
|
testmu(reltol, n, np, p, r,z, props_, condition);
|
|
|
|
testb(reltol,n,np,p,r,z,props_,condition);
|
|
|
|
testrsSat(reltol,n,np,p,props_);
|
|
|
|
testrvSat(reltol,n,np,p,props_);
|
|
|
|
|
|
}
|
|
|
|
BOOST_AUTO_TEST_CASE(test_wetgas)
|
|
{
|
|
|
|
|
|
// read eclipse deck
|
|
const string filename = "wetgas.DATA";
|
|
cout << "Reading deck: " << filename << endl;
|
|
const EclipseGridParser deck (filename);
|
|
|
|
// setup pvt interface
|
|
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
|
|
std::vector<std::shared_ptr<SinglePvtInterface> > 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;
|
|
|
|
// the tolerance for acceptable difference in values
|
|
const double reltol = 1e-9;
|
|
|
|
std::vector<double> p(n);
|
|
std::vector<double> r(n);
|
|
std::vector<PhasePresence> condition(n);
|
|
std::vector<double> z(n * np);
|
|
|
|
|
|
// Used for forward difference calculations
|
|
const double h_p = 1e4;
|
|
const double h_r = 1e-7;
|
|
|
|
// saturated
|
|
p[0] = 1e7;
|
|
p[1] = p[0] + h_p;
|
|
p[2] = p[0];
|
|
|
|
r[0] = 5e-5;
|
|
r[1] = 5e-5;
|
|
r[2] = 5e-5 + h_r;
|
|
|
|
condition[0].setFreeOil();
|
|
condition[1].setFreeOil();
|
|
condition[2].setFreeOil();
|
|
|
|
|
|
// undersaturated
|
|
p[3] = p[0];
|
|
p[4] = p[1];
|
|
p[5] = p[2];
|
|
|
|
r[3] = 1e-5;
|
|
r[4] = 1e-5;
|
|
r[5] = 1e-5 +h_r;
|
|
|
|
// Corresponing z factors, used to compare with the [p,z] interface
|
|
for (int i = 0; i < n; ++i) {
|
|
z[0+i*np] = 0; z[1+i*np] = r[i];
|
|
z[2+i*np] = 1;
|
|
|
|
}
|
|
|
|
testmu(reltol, n, np, p, r,z, props_, condition);
|
|
|
|
testb(reltol,n,np,p,r,z,props_,condition);
|
|
|
|
testrsSat(reltol,n,np,p,props_);
|
|
|
|
testrvSat(reltol,n,np,p,props_);
|
|
|
|
}
|