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

Introduced posibility to change number of sample points for pvt.

Browse Source 
Did change the PVTW calculation so derivatives are exact.
Extended the test functions for pvt and relperm
This commit is contained in:
Halvor Møll Nilsen 2012-08-31 17:01:07 +02:00
parent 7b1501cead
commit ed2aa9da38
11 changed files with 113 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
CMakeLists.txt
View File

@ -3,7 +3,9 @@ project (opm-core)
enable_language(Fortran)
set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR})
SET(CMAKE_BUILD_TYPE "debug")
SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -D_GLIBCXX_DEBUG -pg -Wall -fopenmp -ggdb")
SET(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -pg -Wall -fopenmp -ggdb")
set(Boost_USE_STATIC_LIBS ON)
set(Boost_USE_MULTITHREADED ON)
@ -58,8 +60,8 @@ target_link_libraries(relperm_test
opmcore
)
set_target_properties(opmcore spu_2p PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
set_target_properties(opmcore sim_2p_incomp_reorder PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
set_target_properties(opmcore sim_wateroil PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
set_target_properties(opmcore pvt_test PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
#set_target_properties(opmcore spu_2p PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
#set_target_properties(opmcore sim_2p_incomp_reorder PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
#set_target_properties(opmcore sim_wateroil PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
#set_target_properties(opmcore pvt_test PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
set_target_properties(opmcore relperm_test PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)

3
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
View File

@ -39,7 +39,8 @@ namespace Opm
const parameter::ParameterGroup& param)
{
rock_.init(deck, grid);
pvt_.init(deck);
int samples = param.getDefault("dead_tab_size", 1025);
pvt_.init(deck, samples);
satprops_.init(deck, grid, param);
if (pvt_.numPhases() != satprops_.numPhases()) {

2
opm/core/fluid/SaturationPropsFromDeck.cpp
View File

@ -138,7 +138,7 @@ namespace Opm
}
// Initialize tables.
int tab_size=param.getDefault("tab_size",200);
int tab_size=param.getDefault("tab_size_kr",200);
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_,tab_size);

6
opm/core/fluid/blackoil/BlackoilPvtProperties.cpp
View File

@ -39,7 +39,7 @@ namespace Opm
}
void BlackoilPvtProperties::init(const EclipseGridParser& deck)
void BlackoilPvtProperties::init(const EclipseGridParser& deck,const int samples)
{
typedef std::vector<std::vector<std::vector<double> > > table_t;
// If we need multiple regions, this class and the SinglePvt* classes must change.
@ -78,7 +78,7 @@ namespace Opm
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
if (deck.hasField("PVDO")) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(deck.getPVDO().pvdo_));
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(deck.getPVDO().pvdo_, samples));
} else if (deck.hasField("PVTO")) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(deck.getPVTO().pvto_));
} else if (deck.hasField("PVCDO")) {
@ -90,7 +90,7 @@ namespace Opm
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
if (deck.hasField("PVDG")) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(deck.getPVDG().pvdg_));
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(deck.getPVDG().pvdg_, samples));
} else if (deck.hasField("PVTG")) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(deck.getPVTG().pvtg_));
} else {

2
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp
View File

@ -47,7 +47,7 @@ namespace Opm
BlackoilPvtProperties();
/// Initialize from deck.
void init(const EclipseGridParser& deck);
void init(const EclipseGridParser& deck,const int samples = 16);
/// Number of active phases.
int numPhases() const;

9
opm/core/fluid/blackoil/SinglePvtConstCompr.hpp
View File

@ -106,13 +106,18 @@ namespace Opm
double* output_B,
double* output_dBdp) const
{
B(n, p, 0, output_B);
//B(n, p, 0, output_B);
if (comp_) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_dBdp[i] = -comp_*output_B[i];
//output_dBdp[i] = -comp_*output_B[i];
double x = comp_*(p[i] - ref_press_);
double d= (1.0 + x + 0.5*x*x);
output_B[i] = ref_B_/d;//(1.0 + x + 0.5*x*x);
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);
}
}

4
opm/core/fluid/blackoil/SinglePvtDead.cpp
View File

@ -32,9 +32,8 @@ namespace Opm
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
SinglePvtDead::SinglePvtDead(const table_t& pvd_table)
SinglePvtDead::SinglePvtDead(const table_t& pvd_table,const int samples)
{
const int region_number = 0;
if (pvd_table.size() != 1) {
@ -51,7 +50,6 @@ namespace Opm
B_inv[i] = 1.0 / pvd_table[region_number][1][i];
visc[i] = pvd_table[region_number][2][i];
}
int samples = 1025;
buildUniformMonotoneTable(press, B_inv, samples, one_over_B_);
buildUniformMonotoneTable(press, visc, samples, viscosity_);

3
opm/core/fluid/blackoil/SinglePvtDead.hpp
View File

@ -37,8 +37,7 @@ namespace Opm
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtDead(const table_t& pvd_table);
SinglePvtDead(const table_t& pvd_table,const int samples = 16);
virtual ~SinglePvtDead();
/// Viscosity as a function of p and z.

13
opm/core/utility/UniformTableLinear.hpp
View File

@ -188,9 +188,14 @@ namespace Opm {
UniformTableLinear<T>
::derivative(const double xparam) const
{
// Implements ClosestValue policy.
// Implements derivative consistent
// with ClosestValue policy for function
double value;
if(xparam > xmax_ || xparam < xmin_){
value=0.0;
}else{
double x = std::min(xparam, xmax_);
x = std::max(x, xmin_);
x = std::max(x, xmin_);
// Lookup is easy since we are uniform in x.
double pos = (x - xmin_)/xdelta_;
@ -200,7 +205,9 @@ namespace Opm {
// We are at xmax_
--left;
}
return (y_values_[left + 1] - y_values_[left])/xdelta_;
value = (y_values_[left + 1] - y_values_[left])/xdelta_;
}
return value;
}

92
tests/pvt_test.cpp
View File

@ -30,34 +30,68 @@
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <iterator>
#include <string>
#include <vector>
int main(int argc, char** argv)
{
using namespace std;
// Parameters.
Opm::parameter::ParameterGroup param(argc, argv);
// Parser.
std::string ecl_file = param.get<std::string>("deck_filename");
std::string input_file = param.get<std::string>("input_filename");
std::string matrix_output = param.get<std::string>("matrix_outout");
std::string matrix_output = param.get<std::string>("matrix_output");
std::string rs_output = param.get<std::string>("rs_output");
std::string b_output = param.get<std::string>("b_output");
std::string mu_output = param.get<std::string>("mu_output");
Opm::EclipseGridParser deck(ecl_file);
UnstructuredGrid grid;
grid.number_of_cells = 1;
grid.global_cell = NULL;
Opm::BlackoilPropertiesFromDeck props(deck, grid);
Opm::BlackoilPropertiesFromDeck props(deck, grid, param);
Opm::BlackoilPvtProperties pvt;
pvt.init(deck);
int samples = param.getDefault("dead_tab_size", 1025);
pvt.init(deck, samples);
const int n = 1;
std::fstream inos(input_file.c_str());
if(!inos.good()){
std::cout << "Could not open :" << input_file << std::endl;
exit(3);
}
std::fstream aos(matrix_output.c_str(), std::fstream::out | std::fstream::trunc);
aos << setiosflags(ios::scientific) << setprecision(12);
if(!(aos.good())){
std::cout << "Could not open"<< matrix_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 rsos(rs_output.c_str(), std::fstream::out | std::fstream::trunc);
rsos << setiosflags(ios::scientific) << setprecision(12);
if(!(rsos.good())){
std::cout << "Could not open"<< rs_output << std::endl;
exit(3);
}
std::fstream muos(mu_output.c_str(), std::fstream::out | std::fstream::trunc);
if(!(muos.good())){
std::cout << "Could not open"<< rs_output << std::endl;
exit(3);
}
int np=props.numPhases();
while(inos.good()){
while((inos.good()) && (!inos.eof())){
double p[n];
double z[np*n];
int cells[n] = { 0 };
@ -65,26 +99,36 @@ int main(int argc, char** argv)
for(int i=0; i < np; ++i){
inos >> z[i];
}
double A[np*np*n];
props.matrix(n, p, z, cells, A, 0);
std::fstream aos(matrix_output.c_str());
std::copy(A, A + np*np, std::ostream_iterator<double>(aos, " "));
std::cout << std::endl;
double b[np];
double dbdp[np];
pvt.dBdp(n, p, z, b, dbdp);
std::fstream bos(b_output.c_str());
std::copy(b, b + np, std::ostream_iterator<double>(bos, " "));
std::copy(b, dbdp + np, std::ostream_iterator<double>(bos, " "));
std::cout << std::endl;
double rs[np];
double drs[np];
//pvt.R(n, p, z, rs);
pvt.dRdp(n, p, z, rs,drs);
std::fstream rsos(rs_output.c_str());
std::copy(b, rs + np, std::ostream_iterator<double>(rsos, " "));
std::copy(b, drs + np, std::ostream_iterator<double>(rsos, " "));
std::cout << std::endl;
if(inos.good()){
double A[np*np*n];
double dA[np*np*n];
props.matrix(n, p, z, cells, A, dA);
std::copy(A, A + np*np, std::ostream_iterator<double>(aos, " "));
std::copy(dA, dA + np*np, std::ostream_iterator<double>(aos, " "));
aos << std::endl;
double mu[np];
//double dmu[np];//not implemented
props.viscosity(n, p, z, cells, mu, 0);
std::copy(mu, mu + np, std::ostream_iterator<double>(muos, " "));
//std::copy(dmu, dmu + np, std::ostream_iterator<double>(muos, " "));
aos << std::endl;
double b[np];
double dbdp[np];
pvt.dBdp(n, p, z, b, dbdp);
std::copy(b, b + np, std::ostream_iterator<double>(bos, " "));
std::copy(dbdp, dbdp + np, std::ostream_iterator<double>(bos, " "));
bos << std::endl;
double rs[np];
double drs[np];
//pvt.R(n, p, z, rs);
pvt.dRdp(n, p, z, rs,drs);
std::copy(rs, rs + np, std::ostream_iterator<double>(rsos, " "));
std::copy(drs, drs + np, std::ostream_iterator<double>(rsos, " "));
rsos << std::endl;
}
}
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";

24
tests/relperm_test.cpp
View File

@ -65,20 +65,20 @@ int main(int argc, char** argv)
while((inos.good()) && (!inos.eof())){
double s[np];
for(int i=0; i < np; ++i){
inos >> s[i];
inos >> s[i];
}
if(inos.good()){
double kr[np];
double dkr[np*np];
int cell[1];
cell[0]=1;
props.relperm(1,s, cell, kr, dkr);
std::copy(s, s + np, std::ostream_iterator<double>(kros, " "));
kros << " ";
std::copy(kr, kr + np, std::ostream_iterator<double>(kros, " "));
kros << " ";
std::copy(dkr, dkr + np*np, std::ostream_iterator<double>(kros, " "));
kros << "\n";
double kr[np];
double dkr[np*np];
int cell[1];
cell[0]=1;
props.relperm(1,s, cell, kr, dkr);
std::copy(s, s + np, std::ostream_iterator<double>(kros, " "));
kros << " ";
std::copy(kr, kr + np, std::ostream_iterator<double>(kros, " "));
kros << " ";
std::copy(dkr, dkr + np*np, std::ostream_iterator<double>(kros, " "));
kros << "\n";
}
}
if (param.anyUnused()) {