Merge pull request #254 from totto82/newfluid2

A test for the blackoil fluid based on the [p,r] interface
This commit is contained in:
Bård Skaflestad
2013-06-10 06:50:37 -07:00
4 changed files with 273 additions and 255 deletions

View File

@@ -147,6 +147,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_column_extract.cpp
tests/test_geom2d.cpp
tests/test_param.cpp
tests/test_blackoilfluid.cpp
)
# originally generated with the command:
@@ -154,6 +155,7 @@ list (APPEND TEST_SOURCE_FILES
list (APPEND TEST_DATA_FILES
tests/extratestdata.xml
tests/testdata.xml
tests/testFluid.DATA
)
# originally generated with the command:

View File

@@ -1,255 +0,0 @@
#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>
#include <iostream>
#include <iterator>
#include <vector>
#include <string>
using namespace Opm;
using namespace std;
// The function object divides a Factor with an element
template <class Type>
class MultValue
{
private:
Type Factor; // The value to multiply by
public:
// Constructor initializes the value to multiply by
MultValue ( const Type& _Val ) : Factor ( _Val ) {
}
// The function call for the element to be multiplied
int operator ( ) ( Type& elem ) const
{
return Factor / elem;
}
};
int main () {
// read parameters from command-line
const string filename = "../../opm-core/tests/not-unit/blackoil/SPE9small.DATA";
cout << "Reading deck: " << filename << endl;
const EclipseGridParser deck (filename);
std::string mu_output = "mu_output";
std::string b_output = "b_output";
std::string rbub_output = "rbub_output";
PhaseUsage phase_usage_;
std::vector<std::tr1::shared_ptr<SinglePvtInterface> > props_;
phase_usage_ = phaseUsageFromDeck(deck);
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
int samples = 0;
std::fstream muos(mu_output.c_str(), std::fstream::out | std::fstream::trunc);
if(!(muos.good())){
std::cout << "Could not open"<< mu_output << std::endl;
exit(3);
}
std::fstream bos(b_output.c_str(), std::fstream::out | std::fstream::trunc);
bos << setiosflags(ios::scientific) << setprecision(12);
if(!(bos.good())){
std::cout << "Could not open"<< b_output << std::endl;
exit(3);
}
std::fstream rbubos(rbub_output.c_str(), std::fstream::out | std::fstream::trunc);
rbubos << setiosflags(ios::scientific) << setprecision(12);
if(!(rbubos.good())){
std::cout << "Could not open"<< rbub_output << std::endl;
exit(3);
}
// 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 {
THROW("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 {
THROW("Input is missing PVDG or PVTG\n");
}
}
int n = 6;
int np = 3; //phase_usage_.num_phases;
double p[n];
double r[n];
double z[np*n];
double mu[n];
double dmudp[n];
double dmudr[n];
double mu_new[n];
double dmudp_diff;
double dmudr_diff;
double dmudp_diff_u;
double dmudr_diff_u;
//double rf[n];
double h = 1;
double rh = 1;
p[0] = 10000000;
p[1] = p[0] + h;
p[2] = 10000000;
p[3] = 10000000;
p[4] = p[0] + h;
p[5] = 10000000;
// saturated
r[0] = 200;
r[1] = 200;
r[2] = 200 + rh;
// undersaturated
r[3] = 50;
r[4] = 50;
r[5] = 50 +rh;
for (int i = 0; i < n; ++i) {
z[0+i*np] = 0; z[1+i*np] = 1;
z[2+i*np] = r[i];
}
// test mu
for (int phase = 1; phase < 2; ++phase) {
props_[phase]->mu(n, p, r, mu_new,dmudp,dmudr);
props_[phase]->mu(n, p, z, mu);
dmudp_diff = (mu_new[1]-mu_new[0])/h;
dmudr_diff = (mu_new[2]-mu_new[0])/rh;
dmudp_diff_u = (mu_new[4]-mu_new[3])/h;
dmudr_diff_u = (mu_new[5]-mu_new[3])/rh;
std::copy(mu,mu + n, std::ostream_iterator<double>(muos, " "));
muos << "\n";
std::copy(mu_new,mu_new + n, std::ostream_iterator<double>(muos, " "));
muos << "\n";
std::copy(dmudp,dmudp + n, std::ostream_iterator<double>(muos, " "));
muos << "\n";
muos << dmudp_diff << " " << dmudp_diff_u << "\n";
std::copy(dmudr,dmudr + n, std::ostream_iterator<double>(muos, " "));
muos << "\n";
muos << dmudr_diff << " " << dmudr_diff_u << "\n";
}
// test b
double b[n];
double B[n];
double invB[n];
double dinvBdp[n];
double dBdp[n];
double dbdr[n];
double dbdp[n];
double dbdp_diff;
double dbdr_diff;
double dbdp_diff_u;
double dbdr_diff_u;
for (int phase = 1; phase < 2; ++phase) {
props_[phase]->b(n, p, r, b,dbdp,dbdr);
//props_[phase]->B(n, p, z, B);
props_[phase]->dBdp(n, p, z, B, dBdp);
dbdp_diff = (b[1]-b[0])/h;
dbdr_diff = (b[2]-b[0])/rh;
dbdp_diff_u = (b[4]-b[3])/h;
dbdr_diff_u = (b[5]-b[3])/rh;
for (int i = 0; i < n; ++i){
invB[i] = 1/B[i];
dinvBdp[i] = -1/pow(B[i],2) * dBdp[i];
}
std::copy(b,b + n, std::ostream_iterator<double>(bos, " "));
bos << "\n";
std::copy(invB,invB + n, std::ostream_iterator<double>(bos, " "));
bos << "\n";
std::copy(dinvBdp,dinvBdp + n, std::ostream_iterator<double>(bos, " "));
bos << "\n";
std::copy(dbdp,dbdp + n, std::ostream_iterator<double>(bos, " "));
bos << "\n";
bos << dbdp_diff << " " << dbdp_diff_u << "\n";
std::copy(dbdr,dbdr + n, std::ostream_iterator<double>(bos, " "));
bos << "\n";
bos << dbdr_diff << " " << dbdr_diff_u << "\n";
}
// test rbub
double rbub[n];
double drbubdp[n];
double drbubdp_diff;
for (int phase = 1; phase < 2; ++phase) {
props_[phase] ->rbub(n,p,rbub,drbubdp);
drbubdp_diff = (rbub[1]-rbub[0])/h;
std::copy(rbub,rbub + n, std::ostream_iterator<double>(rbubos, " "));
rbubos << "\n";
std::copy(drbubdp,drbubdp + n, std::ostream_iterator<double>(rbubos, " "));
rbubos << drbubdp_diff;
rbubos << "\n";
}
}

44
tests/testFluid.DATA Normal file
View File

@@ -0,0 +1,44 @@
OIL
WATER
GAS
FIELD
PVTO
-- Rs Pbub Bo Vo
.0 14.7 1.0000 1.20 /
.165 400. 1.0120 1.17 /
.335 800. 1.0255 1.14 /
.500 1200. 1.0380 1.11 /
.665 1600. 1.0510 1.08 /
.828 2000. 1.0630 1.06 /
.985 2400. 1.0750 1.03 /
1.130 2800. 1.0870 1.00 /
1.270 3200. 1.0985 .98 /
1.390 3600. 1.1100 .95 /
1.500 4000. 1.1200 .94
5000. 1.1189 .94 /
/
PVDG
-- Pg Bg Vg
14.7 178.08 .0125
400. 5.4777 .0130
800. 2.7392 .0135
1200. 1.8198 .0140
1600. 1.3648 .0145
2000. 1.0957 .0150
2400. 0.9099 .0155
2800. 0.7799 .0160
3200. 0.6871 .0165
3600. 0.6035 .0170
4000. 0.5432 .0175 /
PVTW
--Depth Bw Comp Vw Cv
3600. 1.0034 1.0E-6 0.96 0.0 /
DENSITY
-- Oil Water Gas
44.98 63.01 0.0702 /

View File

@@ -0,0 +1,227 @@
#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 <iostream>
#include <iterator>
#include <vector>
#include <string>
using namespace Opm;
using namespace std;
BOOST_AUTO_TEST_CASE(test_blackoilfluid)
{
// read eclipse deck
const string filename = "testFluid.DATA";
cout << "Reading deck: " << filename << endl;
const EclipseGridParser deck (filename);
// setup pvt interface
std::vector<std::tr1::shared_ptr<SinglePvtInterface> > props_;
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
int samples = 0;
// 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 {
THROW("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 {
THROW("Input is missing PVDG or PVTG\n");
}
}
// 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<double> z(n * np);
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;
// Used for forward difference calculations
const double h_p = 100000;
const double h_r = 1;
// saturated
p[0] = 10000000;
p[1] = p[0] + h_p;
p[2] = p[0];
r[0] = 200;
r[1] = 200;
r[2] = 200 + h_r;
// 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];
}
// test mu
for (int phase = 1; phase < 2; ++phase) {
props_[phase]->mu(n, &p[0], &r[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])/h_p;
dmudr_diff = (mu_new[2]-mu_new[0])/h_r;
dmudp_diff_u = (mu_new[4]-mu_new[3])/h_p;
dmudr_diff_u = (mu_new[5]-mu_new[3])/h_r;
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);
}
// 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 = 1; phase < 2; ++phase) {
props_[phase]->b(n, &p[0], &r[0], &b[0], &dbdp[0], &dbdr[0]);
//props_[phase]->B(n, p, z, B);
props_[phase]->dBdp(n, &p[0], &z[0], &B[0], &dBdp[0]);
dbdp_diff = (b[1]-b[0])/h_p;
dbdr_diff = (b[2]-b[0])/h_r;
dbdp_diff_u = (b[4]-b[3])/h_p;
dbdr_diff_u = (b[5]-b[3])/h_r;
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);
}
// test bublepoint pressure
std::vector<double> rbub(n);
std::vector<double> drbubdp(n);
double drbubdp_diff;
double drbubdp_diff_u;
for (int phase = 1; phase < 2; ++phase) {
props_[phase] ->rbub(n, &p[0], &rbub[0], &drbubdp[0]);
drbubdp_diff = (rbub[1]-rbub[0])/h_p;
drbubdp_diff_u = (rbub[4]-rbub[3])/h_p;
// saturated case
BOOST_CHECK_CLOSE(drbubdp_diff,drbubdp[0], reltol);
// unsaturad case
BOOST_CHECK_CLOSE(drbubdp_diff_u,drbubdp[3], reltol);
}
}