diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index 8807650a6..4ca597a4d 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -20,14 +20,15 @@ #include "config.h" #include #include +#include #include #include +#include #include #include namespace Opm { - BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState, const UnstructuredGrid& grid, @@ -133,16 +134,15 @@ namespace Opm if (init_rock){ rock_.init(eclState, number_of_cells, global_cell, cart_dims); } - pvt_.init(deck, eclState, /*numSamples=*/0); + phaseUsage_ = phaseUsageFromDeck(deck); + initSurfaceDensities_(deck); + oilPvt_.initFromDeck(deck, eclState); + gasPvt_.initFromDeck(deck, eclState); + waterPvt_.initFromDeck(deck, eclState); SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); ptr->init(phaseUsageFromDeck(deck), materialLawManager); satprops_.reset(ptr); - - if (pvt_.numPhases() != satprops_->numPhases()) { - OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); - } } inline void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr deck, @@ -162,8 +162,11 @@ namespace Opm rock_.init(eclState, number_of_cells, global_cell, cart_dims); } - const int pvt_samples = param.getDefault("pvt_tab_size", -1); - pvt_.init(deck, eclState, pvt_samples); + phaseUsage_ = phaseUsageFromDeck(deck); + initSurfaceDensities_(deck); + oilPvt_.initFromDeck(deck, eclState); + gasPvt_.initFromDeck(deck, eclState); + waterPvt_.initFromDeck(deck, eclState); // Unfortunate lack of pointer smartness here... std::string threephase_model = param.getDefault("threephase_model", "gwseg"); @@ -175,18 +178,12 @@ namespace Opm = new SaturationPropsFromDeck(); ptr->init(phaseUsageFromDeck(deck), materialLawManager); satprops_.reset(ptr); - - if (pvt_.numPhases() != satprops_->numPhases()) { - OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); - } } BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck() { } - /// \return D, the number of spatial dimensions. int BlackoilPropertiesFromDeck::numDimensions() const { @@ -219,13 +216,13 @@ namespace Opm /// \return P, the number of phases (also the number of components). int BlackoilPropertiesFromDeck::numPhases() const { - return pvt_.numPhases(); + return phaseUsage_.num_phases; } /// \return Object describing the active phases. PhaseUsage BlackoilPropertiesFromDeck::phaseUsage() const { - return pvt_.phaseUsage(); + return phaseUsage_; } /// \param[in] n Number of data points. @@ -244,16 +241,45 @@ namespace Opm double* mu, double* dmudp) const { - if (dmudp) { - OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::viscosity() -- derivatives of viscosity not yet implemented."); - } else { - const int *cellPvtTableIdx = cellPvtRegionIndex(); - assert(cellPvtTableIdx != 0); - std::vector pvtTableIdx(n); - for (int i = 0; i < n; ++ i) - pvtTableIdx[i] = cellPvtTableIdx[cells[i]]; + const auto& pu = phaseUsage(); - pvt_.mu(n, &pvtTableIdx[0], p, T, z, mu); + enum PressureEvalTag {}; + typedef Opm::LocalAd::Evaluation LadEval; + + LadEval pLad = 0.0; + LadEval TLad = 0.0; + LadEval RsLad = 0.0; + LadEval RvLad = 0.0; + LadEval muLad = 0.0; + + pLad.derivatives[0] = 1.0; + + for (int i = 0; i < n; ++ i) { + int cellIdx = cells[i]; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + pLad.value = p[i]; + TLad.value = T[i]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + muLad = waterPvt_.viscosity(pvtRegionIdx, TLad, pLad); + int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua]; + mu[offset] = muLad.value; + dmudp[offset] = muLad.derivatives[0]; + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + muLad = oilPvt_.viscosity(pvtRegionIdx, TLad, pLad, RsLad); + int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]; + mu[offset] = muLad.value; + dmudp[offset] = muLad.derivatives[0]; + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + muLad = gasPvt_.viscosity(pvtRegionIdx, TLad, pLad, RvLad); + int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]; + mu[offset] = muLad.value; + dmudp[offset] = muLad.derivatives[0]; + } } } @@ -278,27 +304,23 @@ namespace Opm { const int np = numPhases(); - const int *cellPvtTableIdx = cellPvtRegionIndex(); - std::vector 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, &pvtTableIdx[0], p, T, z, &B_[0], &dB_[0]); - pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]); + + this->compute_dBdp_(n, p, T, z, cells, &B_[0], &dB_[0]); + this->compute_dRdp_(n, p, T, z, cells, &R_[0], &dR_[0]); } else { - pvt_.B(n, &pvtTableIdx[0], p, T, z, &B_[0]); - pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]); + this->compute_B_(n, p, T, z, cells, &B_[0]); + this->compute_R_(n, p, T, z, cells, &R_[0]); } - const int* phase_pos = pvt_.phasePosition(); - bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] && - pvt_.phaseUsed()[BlackoilPhases::Vapour]; - const int o = phase_pos[BlackoilPhases::Liquid]; - const int g = phase_pos[BlackoilPhases::Vapour]; + const auto& pu = phaseUsage(); + bool oil_and_gas = pu.phase_pos[BlackoilPhases::Liquid] && + pu.phase_pos[BlackoilPhases::Vapour]; + const int o = pu.phase_pos[BlackoilPhases::Liquid]; + const int g = pu.phase_pos[BlackoilPhases::Vapour]; // Compute A matrix // #pragma omp parallel for @@ -360,6 +382,276 @@ namespace Opm } } + void BlackoilPropertiesFromDeck::compute_B_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* B) const + { + const auto& pu = phaseUsage(); + + typedef double LadEval; + + LadEval pLad = 0.0; + LadEval TLad = 0.0; + LadEval RsLad = 0.0; + LadEval RvLad = 0.0; + + for (int i = 0; i < n; ++ i) { + int cellIdx = cells[i]; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + pLad = p[i]; + TLad = T[i]; + + int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid]; + int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour]; + int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + LadEval BLad = 1.0/waterPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + + B[waterOffset] = BLad; + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + double currentRs = 0.0; + double maxRs = 0.0; + if (pu.phase_used[BlackoilPhases::Vapour]) { + currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset]; + maxRs = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad); + } + LadEval BLad; + if (currentRs >= maxRs) { + BLad = 1.0/oilPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + } + else { + RsLad = currentRs; + BLad = 1.0/oilPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad); + } + + B[oilOffset] = BLad; + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + double currentRv = 0.0; + double maxRv = 0.0; + if (pu.phase_used[BlackoilPhases::Liquid]) { + currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset]; + maxRv = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad); + } + LadEval BLad; + if (currentRv >= maxRv) { + BLad = 1.0/gasPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + } + else { + RvLad = currentRv; + BLad = 1.0/gasPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RvLad); + } + + B[gasOffset] = BLad; + } + } + } + + void BlackoilPropertiesFromDeck::compute_dBdp_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* B, + double* dBdp) const + { + const auto& pu = phaseUsage(); + + enum PressureEvalTag {}; + typedef Opm::LocalAd::Evaluation LadEval; + + LadEval pLad = 0.0; + LadEval TLad = 0.0; + LadEval RsLad = 0.0; + LadEval RvLad = 0.0; + + pLad.derivatives[0] = 1.0; + + for (int i = 0; i < n; ++ i) { + int cellIdx = cells[i]; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + pLad.value = p[i]; + TLad.value = T[i]; + + int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid]; + int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour]; + int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + LadEval BLad = 1.0/waterPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + + B[waterOffset] = BLad.value; + dBdp[waterOffset] = BLad.derivatives[0]; + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + double currentRs = 0.0; + double maxRs = 0.0; + if (pu.phase_used[BlackoilPhases::Vapour]) { + currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset]; + maxRs = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad.value, pLad.value); + } + LadEval BLad; + if (currentRs >= maxRs) { + BLad = 1.0/oilPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + } + else { + RsLad.value = currentRs; + BLad = 1.0/oilPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad); + } + + B[oilOffset] = BLad.value; + dBdp[oilOffset] = BLad.derivatives[0]; + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + double currentRv = 0.0; + double maxRv = 0.0; + if (pu.phase_used[BlackoilPhases::Liquid]) { + currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset]; + maxRv = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad.value, pLad.value); + } + LadEval BLad; + if (currentRv >= maxRv) { + BLad = 1.0/gasPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + } + else { + RvLad.value = currentRv; + BLad = 1.0/gasPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RvLad); + } + + B[gasOffset] = BLad.value; + dBdp[gasOffset] = BLad.derivatives[0]; + } + } + } + + void BlackoilPropertiesFromDeck::compute_R_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* R) const + { + const auto& pu = phaseUsage(); + + typedef double LadEval; + + LadEval pLad = 0.0; + LadEval TLad = 0.0; + + for (int i = 0; i < n; ++ i) { + int cellIdx = cells[i]; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + pLad = p[i]; + TLad = T[i]; + + int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid]; + int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour]; + int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + R[waterOffset] = 0.0; // water is always immiscible! + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + LadEval RsSatLad = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad); + + double currentRs = 0.0; + if (pu.phase_used[BlackoilPhases::Vapour]) { + currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset]; + } + + RsSatLad = std::min(RsSatLad, currentRs); + + R[oilOffset] = RsSatLad; + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + LadEval RvSatLad = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad); + + double currentRv = 0.0; + if (pu.phase_used[BlackoilPhases::Liquid]) { + currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset]; + } + + RvSatLad = std::min(RvSatLad, currentRv); + + R[gasOffset] = RvSatLad; + } + } + } + + void BlackoilPropertiesFromDeck::compute_dRdp_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* R, + double* dRdp) const + { + const auto& pu = phaseUsage(); + + enum PressureEvalTag {}; + typedef Opm::LocalAd::Evaluation LadEval; + typedef Opm::MathToolbox Toolbox; + + LadEval pLad = 0.0; + LadEval TLad = 0.0; + + pLad.derivatives[0] = 1.0; + + for (int i = 0; i < n; ++ i) { + int cellIdx = cells[i]; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + pLad.value = p[i]; + TLad.value = T[i]; + + int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid]; + int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour]; + int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + R[waterOffset] = 0.0; // water is always immiscible! + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + LadEval RsSatLad = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad); + + LadEval currentRs = 0.0; + if (pu.phase_used[BlackoilPhases::Vapour]) { + currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset]; + } + + RsSatLad = Toolbox::min(RsSatLad, currentRs); + + R[oilOffset] = RsSatLad.value; + dRdp[oilOffset] = RsSatLad.derivatives[0]; + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + LadEval RvSatLad = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad); + + LadEval currentRv = 0.0; + if (pu.phase_used[BlackoilPhases::Liquid]) { + currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset]; + } + + RvSatLad = Toolbox::min(RvSatLad, currentRv); + + R[gasOffset] = RvSatLad.value; + dRdp[gasOffset] = RvSatLad.derivatives[0]; + } + } + } + /// \param[in] n Number of data points. /// \param[in] A Array of nP^2 values, where the P^2 values for a cell give the /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices @@ -376,8 +668,7 @@ namespace Opm // #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); + const double *sdens = surfaceDensity(cellIdx); for (int phase = 0; phase < np; ++phase) { rho[np*i + phase] = 0.0; for (int comp = 0; comp < np; ++comp) { @@ -391,8 +682,37 @@ namespace Opm /// \return Array of P density values. const double* BlackoilPropertiesFromDeck::surfaceDensity(int cellIdx) const { + const auto& pu = phaseUsage(); int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx); - return pvt_.surfaceDensities(pvtRegionIdx); + return &surfaceDensities_[pvtRegionIdx*pu.num_phases]; + } + + void BlackoilPropertiesFromDeck::initSurfaceDensities_(Opm::DeckConstPtr deck) + { + const auto& pu = phaseUsage(); + int np = pu.num_phases; + int numPvtRegions = 1; + if (deck->hasKeyword("TABDIMS")) { + const auto& tabdimsKeyword = deck->getKeyword("TABDIMS"); + numPvtRegions = tabdimsKeyword.getRecord(0).getItem("NTPVT").template get(0); + } + + const auto& densityKeyword = deck->getKeyword("DENSITY"); + + surfaceDensities_.resize(np*numPvtRegions); + for (int pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++pvtRegionIdx) { + if (pu.phase_used[BlackoilPhases::Aqua]) + surfaceDensities_[np*pvtRegionIdx + pu.phase_pos[BlackoilPhases::Aqua]] = + densityKeyword.getRecord(pvtRegionIdx).getItem("WATER").getSIDouble(0); + + if (pu.phase_used[BlackoilPhases::Liquid]) + surfaceDensities_[np*pvtRegionIdx + pu.phase_pos[BlackoilPhases::Liquid]] = + densityKeyword.getRecord(pvtRegionIdx).getItem("OIL").getSIDouble(0); + + if (pu.phase_used[BlackoilPhases::Vapour]) + surfaceDensities_[np*pvtRegionIdx + pu.phase_pos[BlackoilPhases::Vapour]] = + densityKeyword.getRecord(pvtRegionIdx).getItem("GAS").getSIDouble(0); + } } /// \param[in] n Number of data points. diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index fb37a7f89..e5c183888 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -23,9 +23,11 @@ #include #include -#include #include #include +#include +#include +#include #include @@ -233,10 +235,19 @@ namespace Opm const double pcow, double & swat); - // return a reference to the "raw" PVT fluid object for a phase. - const PvtInterface& pvt(int phaseIdx) const + const OilPvtMultiplexer& oilPvt() const { - return pvt_.pvt(phaseIdx); + return oilPvt_; + } + + const GasPvtMultiplexer& gasPvt() const + { + return gasPvt_; + } + + const WaterPvtMultiplexer& waterPvt() const + { + return waterPvt_; } private: @@ -247,6 +258,38 @@ namespace Opm return pvtTableIdx[cellIdx]; } + void initSurfaceDensities_(Opm::DeckConstPtr deck); + + void compute_B_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* B) const; + + void compute_dBdp_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* B, + double* dBdp) const; + + void compute_R_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* R) const; + + void compute_dRdp_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* R, + double* dRdp) const; + void init(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState, std::shared_ptr materialLawManager, @@ -265,10 +308,14 @@ namespace Opm bool init_rock); RockFromDeck rock_; + PhaseUsage phaseUsage_; std::vector cellPvtRegionIdx_; - BlackoilPvtProperties pvt_; + OilPvtMultiplexer oilPvt_; + GasPvtMultiplexer gasPvt_; + WaterPvtMultiplexer waterPvt_; std::shared_ptr materialLawManager_; std::shared_ptr satprops_; + std::vector surfaceDensities_; mutable std::vector B_; mutable std::vector dB_; mutable std::vector R_; diff --git a/opm/core/props/IncompPropertiesFromDeck.hpp b/opm/core/props/IncompPropertiesFromDeck.hpp index e71ecb781..4932ed805 100644 --- a/opm/core/props/IncompPropertiesFromDeck.hpp +++ b/opm/core/props/IncompPropertiesFromDeck.hpp @@ -23,8 +23,8 @@ #include #include #include -#include #include +#include #include struct UnstructuredGrid; diff --git a/tests/norne_pvt.data b/tests/norne_pvt.data index 10de78063..bdb7af927 100644 --- a/tests/norne_pvt.data +++ b/tests/norne_pvt.data @@ -5,11 +5,22 @@ -- Copyright (C) 2015 Statoil +RUNSPEC TABDIMS --ntsfun ntpvt nssfun nppvt ntfip nrpvt ntendp 2 2 33 60 16 60 / +DIMENS +1 1 1 / + +GRID + +PROPS + +DENSITY + 859.5 1033.0 0.854 / Justert 22/7 + 860.04 1033.0 0.853 / Justert 22/7 PVTG diff --git a/tests/test_norne_pvt.cpp b/tests/test_norne_pvt.cpp index 39a9c868a..57e28fcd9 100644 --- a/tests/test_norne_pvt.cpp +++ b/tests/test_norne_pvt.cpp @@ -39,7 +39,7 @@ #include #include -#include +#include using namespace Opm; @@ -57,19 +57,9 @@ using namespace Opm; further semantic meaning. */ - -void check_vectors( const std::vector& v1 , const std::vector& v2) { - double tol = 1e-5; - for (decltype(v1.size()) i = 0; i < v1.size(); i++) { - BOOST_CHECK_CLOSE( v1[i] , v2[i] , tol ); - } -} - - - -void verify_norne_oil_pvt_region2(const TableManager& tableManager) { - auto pvtoTables = tableManager.getPvtoTables(); - PvtLiveOil oilPvt(pvtoTables); +void verify_norne_oil_pvt_region1(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState) { + Opm::LiveOilPvt oilPvt; + oilPvt.initFromDeck(deck, eclState); std::vector rs = {33, 33, 43, 43, @@ -113,33 +103,37 @@ void verify_norne_oil_pvt_region2(const TableManager& tableManager) { { std::vector tableIndex(P.size() , 0); - std::vector mu(P.size()); - std::vector dmudp(P.size()); - std::vector dmudr(P.size()); - - std::vector b(P.size()); - std::vector dbdp(P.size()); - std::vector dbdr(P.size()); - - + // convert the pressures to SI units (bar to Pascal) for (auto& value : P) value *= Metric::Pressure; + // convert the gas dissolution factors to SI units for (auto& value : rs) value *= Metric::GasDissolutionFactor; - oilPvt.mu( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , mu.data() , dmudp.data() , dmudr.data()); - oilPvt.b( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , b.data() , dbdp.data() , dbdr.data()); + for (unsigned i = 0; i < P.size(); ++i) { + double mu; + double b; + double RsSat = oilPvt.saturatedGasDissolutionFactor(/*tableIndex=*/0, /*T=*/273.15, P[i]); + if (rs[i] >= RsSat) { + mu = oilPvt.saturatedViscosity(/*tableIndex=*/0, /*T=*/273.15, P[i]); + b = oilPvt.saturatedInverseFormationVolumeFactor(/*tableIndex=*/0, /*T=*/273.15, P[i]); + } + else { + mu = oilPvt.viscosity(/*tableIndex=*/0, /*T=*/273.15, P[i], rs[i]); + b = oilPvt.inverseFormationVolumeFactor(/*tableIndex=*/0, /*T=*/273.15, P[i], rs[i]); + } - check_vectors( mu , mu_expected ); - check_vectors( b , b_expected ); + BOOST_CHECK_CLOSE( mu , mu_expected[i], 1e-5 ); + BOOST_CHECK_CLOSE( b , b_expected[i], 1e-5 ); + } } } -void verify_norne_oil_pvt_region1(const TableManager& tableManager) { - auto pvtoTables = tableManager.getPvtoTables(); - PvtLiveOil oilPvt(pvtoTables); +void verify_norne_oil_pvt_region2(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState) { + Opm::LiveOilPvt oilPvt; + oilPvt.initFromDeck(deck, eclState); std::vector rs = {21 , 21, 30 , 30, @@ -254,50 +248,41 @@ void verify_norne_oil_pvt_region1(const TableManager& tableManager) { 0.57289091037, 0.56019050084, 0.55474601877, 0.55809201119, 0.54526832277}; - { - std::vector tableIndex(P.size() , 1); + // convert the pressures to SI units (bar to Pascal) + for (auto& value : P) + value *= Metric::Pressure; - std::vector mu(P.size()); - std::vector dmudp(P.size()); - std::vector dmudr(P.size()); + // convert the gas dissolution factors to SI units + for (auto& value : rs) + value *= Metric::GasDissolutionFactor; - std::vector b(P.size()); - std::vector dbdp(P.size()); - std::vector dbdr(P.size()); + for (unsigned i = 0; i < P.size(); ++i) { + double mu; + double b; + double RsSat = oilPvt.saturatedGasDissolutionFactor(/*tableIndex=*/1, /*T=*/273.15, P[i]); + if (rs[i] >= RsSat) { + mu = oilPvt.saturatedViscosity(/*tableIndex=*/1, /*T=*/273.15, P[i]); + b = oilPvt.saturatedInverseFormationVolumeFactor(/*tableIndex=*/1, /*T=*/273.15, P[i]); + } + else { + mu = oilPvt.viscosity(/*tableIndex=*/1, /*T=*/273.15, P[i], rs[i]); + b = oilPvt.inverseFormationVolumeFactor(/*tableIndex=*/1, /*T=*/273.15, P[i], rs[i]); + } - - for (auto& value : P) - value *= Metric::Pressure; - - for (auto& value : rs) - value *= Metric::GasDissolutionFactor; - - oilPvt.mu( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , mu.data() , dmudp.data() , dmudr.data()); - oilPvt.b( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , b.data() , dbdp.data() , dbdr.data()); - - check_vectors( mu , mu_expected ); - check_vectors( b , b_expected ); + BOOST_CHECK_CLOSE( mu , mu_expected[i], 1e-5 ); + BOOST_CHECK_CLOSE( b , b_expected[i], 1e-5 ); } } - - - - -TableManager loadTables( const std::string& deck_file) { +BOOST_AUTO_TEST_CASE( Test_Norne_PVT) { Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }}); Opm::ParserPtr parser(new Parser()); + std::shared_ptr deck; + deck = parser->parseFile("norne_pvt.data", parseMode); - deck = parser->parseFile(deck_file, parseMode); - return TableManager(*deck); -} - - - -BOOST_AUTO_TEST_CASE( Test_Norne_PVT) { - TableManager tableManager = loadTables( "norne_pvt.data" ); - - verify_norne_oil_pvt_region1( tableManager ); - verify_norne_oil_pvt_region2( tableManager ); + Opm::EclipseStateConstPtr eclState(new EclipseState(deck, parseMode)); + + verify_norne_oil_pvt_region1( deck, eclState ); + verify_norne_oil_pvt_region2( deck, eclState ); }