From cbb5910a9323ef9e0e03a3ae262f6b9f38066a79 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 13 May 2014 12:59:11 +0200 Subject: [PATCH] incorperate the review comments/decisions for multi-region PVT the largest change is that all classes below opm/core/props/pvt take the PVT region index as an argument, the higher-level ones (i.e., BlackoilProps*) take cell indices. --- opm/core/pressure/CompressibleTpfa.cpp | 4 +- opm/core/props/BlackoilPropertiesBasic.cpp | 4 +- opm/core/props/BlackoilPropertiesBasic.hpp | 10 ++- opm/core/props/BlackoilPropertiesFromDeck.cpp | 36 ++++++--- opm/core/props/BlackoilPropertiesFromDeck.hpp | 17 +++- .../props/BlackoilPropertiesFromDeck_impl.hpp | 11 ++- .../props/BlackoilPropertiesInterface.hpp | 8 +- opm/core/props/pvt/BlackoilPvtProperties.cpp | 43 +++++----- opm/core/props/pvt/BlackoilPvtProperties.hpp | 38 ++++----- opm/core/props/pvt/PvtConstCompr.hpp | 64 +++++++++------ opm/core/props/pvt/PvtDead.cpp | 37 +++++---- opm/core/props/pvt/PvtDead.hpp | 44 ++++++++--- opm/core/props/pvt/PvtDeadSpline.cpp | 33 ++++---- opm/core/props/pvt/PvtDeadSpline.hpp | 24 +++++- opm/core/props/pvt/PvtInterface.hpp | 43 +++++++--- opm/core/props/pvt/PvtLiveGas.cpp | 42 ++++++---- opm/core/props/pvt/PvtLiveGas.hpp | 24 +++++- opm/core/props/pvt/PvtLiveOil.cpp | 78 ++++++++++--------- opm/core/props/pvt/PvtLiveOil.hpp | 24 +++++- opm/core/props/pvt/PvtPropertiesBasic.hpp | 2 + .../props/pvt/PvtPropertiesIncompFromDeck.cpp | 4 +- opm/core/simulator/EquilibrationHelpers.hpp | 2 +- opm/core/simulator/initStateEquil.hpp | 2 - opm/core/simulator/initState_impl.hpp | 2 +- ...sportSolverCompressibleTwophaseReorder.cpp | 2 +- opm/core/utility/miscUtilitiesBlackoil.cpp | 3 +- tests/test_blackoilfluid.cpp | 64 +++++++-------- 27 files changed, 417 insertions(+), 248 deletions(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 606e04ba..24ac6438 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -238,7 +238,7 @@ namespace Opm const int cell = wells_->well_cells[j]; const double cell_depth = grid_.cell_centroids[dim * cell + dim - 1]; props_.matrix(1, &state.pressure()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0); - props_.density(1, &A[0], &rho[0]); + props_.density(1, &A[0], &cell, &rho[0]); for (int phase = 0; phase < np; ++phase) { const double s_phase = state.saturation()[np*cell + phase]; wellperf_wdp_[j] += s_phase*rho[phase]*grav*(cell_depth - ref_depth); @@ -380,7 +380,7 @@ namespace Opm // Gravity contribution, gravcontrib = rho*(face_z - cell_z) [per phase]. if (grav != 0.0) { const double depth_diff = face_depth - grid_.cell_centroids[c[j]*dim + dim - 1]; - props_.density(1, &cell_A_[np*np*c[j]], &gravcontrib[j][0]); + props_.density(1, &cell_A_[np*np*c[j]], &c[j], &gravcontrib[j][0]); for (int p = 0; p < np; ++p) { gravcontrib[j][p] *= depth_diff*grav; } diff --git a/opm/core/props/BlackoilPropertiesBasic.cpp b/opm/core/props/BlackoilPropertiesBasic.cpp index 8700af11..ab9b238b 100644 --- a/opm/core/props/BlackoilPropertiesBasic.cpp +++ b/opm/core/props/BlackoilPropertiesBasic.cpp @@ -157,9 +157,11 @@ namespace Opm /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices /// are assumed to be in Fortran order, and are typically the result /// of a call to the method matrix(). + /// \param[in] cells The index of the grid cell of each data point. /// \param[out] rho Array of nP density values, array must be valid before calling. void BlackoilPropertiesBasic::density(const int n, const double* A, + const int* /*cells*/, double* rho) const { const int np = numPhases(); @@ -177,7 +179,7 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of P density values. - const double* BlackoilPropertiesBasic::surfaceDensity() const + const double* BlackoilPropertiesBasic::surfaceDensity(int /*cellIdx*/) const { return pvt_.surfaceDensities(); } diff --git a/opm/core/props/BlackoilPropertiesBasic.hpp b/opm/core/props/BlackoilPropertiesBasic.hpp index ad84b1a3..b85991ed 100644 --- a/opm/core/props/BlackoilPropertiesBasic.hpp +++ b/opm/core/props/BlackoilPropertiesBasic.hpp @@ -59,6 +59,11 @@ namespace Opm /// \return N, the number of cells. virtual int numCells() const; + /// Return an array containing the PVT table index for each + /// grid cell + virtual const int* cellPvtRegionIndex() const + { return NULL; } + /// \return Array of N porosity values. virtual const double* porosity() const; @@ -114,14 +119,17 @@ namespace Opm /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices /// are assumed to be in Fortran order, and are typically the result /// of a call to the method matrix(). + /// \param[in] cells The index of the grid cell of each data point. /// \param[out] rho Array of nP density values, array must be valid before calling. virtual void density(const int n, const double* A, + const int* cells, double* rho) const; /// Densities of stock components at surface conditions. + /// \param[in] cellIdx The index of the cell for which the surface density ought to be calculated /// \return Array of P density values. - virtual const double* surfaceDensity() const; + virtual const double* surfaceDensity(int cellIdx = 0) const; /// \param[in] n Number of data points. /// \param[in] s Array of nP saturation values. diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index 5a781031..71246f9c 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -96,14 +96,19 @@ namespace Opm void BlackoilPropertiesFromDeck::viscosity(const int n, const double* p, const double* z, - const int* /*cells*/, + const int* cells, double* mu, double* dmudp) const { if (dmudp) { OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::viscosity() -- derivatives of viscosity not yet implemented."); } else { - pvt_.mu(n, p, z, mu); + const int *cellPvtTableIdx = cellPvtRegionIndex(); + std::vector pvtTableIdx(n); + for (int i = 0; i < n; ++ i) + pvtTableIdx[i] = cellPvtTableIdx[cells[i]]; + + pvt_.mu(n, &pvtTableIdx[0], p, z, mu); } } @@ -120,21 +125,27 @@ namespace Opm void BlackoilPropertiesFromDeck::matrix(const int n, const double* p, const double* z, - const int* /*cells*/, + const int* cells, double* A, double* dAdp) const { 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, p, z, &B_[0], &dB_[0]); - pvt_.dRdp(n, p, z, &R_[0], &dR_[0]); + pvt_.dBdp(n, &pvtTableIdx[0], p, z, &B_[0], &dB_[0]); + pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]); } else { - pvt_.B(n, p, z, &B_[0]); - pvt_.R(n, p, z, &R_[0]); + pvt_.B(n, &pvtTableIdx[0], p, z, &B_[0]); + pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]); } const int* phase_pos = pvt_.phasePosition(); bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] && @@ -207,15 +218,19 @@ namespace Opm /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices /// are assumed to be in Fortran order, and are typically the result /// of a call to the method matrix(). + /// \param[in] cells The index of the grid cell of each data point. /// \param[out] rho Array of nP density values, array must be valid before calling. void BlackoilPropertiesFromDeck::density(const int n, const double* A, + const int* cells, double* rho) const { const int np = numPhases(); - const double* sdens = pvt_.surfaceDensities(); // #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); for (int phase = 0; phase < np; ++phase) { rho[np*i + phase] = 0.0; for (int comp = 0; comp < np; ++comp) { @@ -227,9 +242,10 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of P density values. - const double* BlackoilPropertiesFromDeck::surfaceDensity() const + const double* BlackoilPropertiesFromDeck::surfaceDensity(int cellIdx) const { - return pvt_.surfaceDensities(); + int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx); + return pvt_.surfaceDensities(pvtRegionIdx); } /// \param[in] n Number of data points. diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index e2919424..448c5cb5 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -97,6 +97,11 @@ namespace Opm /// \return N, the number of cells. virtual int numCells() const; + /// Return an array containing the PVT table index for each + /// grid cell + virtual const int* cellPvtRegionIndex() const + { return &cellPvtRegionIdx_[0]; } + /// \return Array of N porosity values. virtual const double* porosity() const; @@ -152,14 +157,16 @@ namespace Opm /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices /// are assumed to be in Fortran order, and are typically the result /// of a call to the method matrix(). + /// \param[in] cells The index of the grid cell of each data point. /// \param[out] rho Array of nP density values, array must be valid before calling. virtual void density(const int n, const double* A, + const int* cells, double* rho) const; /// Densities of stock components at surface conditions. /// \return Array of P density values. - virtual const double* surfaceDensity() const; + virtual const double* surfaceDensity(int cellIdx = 0) const; /// \param[in] n Number of data points. /// \param[in] s Array of nP saturation values. @@ -206,6 +213,13 @@ namespace Opm double* smax) const; private: + int getTableIndex_(const int* pvtTableIdx, int cellIdx) const + { + if (!pvtTableIdx) + return 0; + return pvtTableIdx[cellIdx]; + } + template void init(Opm::DeckConstPtr deck, int number_of_cells, @@ -224,6 +238,7 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock); RockFromDeck rock_; + std::vector cellPvtRegionIdx_; BlackoilPvtProperties pvt_; std::unique_ptr satprops_; mutable std::vector B_; diff --git a/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp b/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp index 713ede5e..586021c4 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp @@ -42,10 +42,14 @@ namespace Opm int dimension, bool init_rock) { + // retrieve the cell specific PVT table index from the deck + // and using the grid... + extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell); + if (init_rock){ rock_.init(deck, number_of_cells, global_cell, cart_dims); } - pvt_.init(deck, /*numSamples=*/0, number_of_cells, global_cell); + pvt_.init(deck, /*numSamples=*/0); SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); satprops_.reset(ptr); @@ -68,13 +72,16 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock) { + // retrieve the cell specific PVT table index from the deck + // and using the grid... + extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell); if(init_rock){ rock_.init(deck, number_of_cells, global_cell, cart_dims); } const int pvt_samples = param.getDefault("pvt_tab_size", 200); - pvt_.init(deck, pvt_samples, number_of_cells, global_cell); + pvt_.init(deck, pvt_samples); // Unfortunate lack of pointer smartness here... const int sat_samples = param.getDefault("sat_tab_size", 200); diff --git a/opm/core/props/BlackoilPropertiesInterface.hpp b/opm/core/props/BlackoilPropertiesInterface.hpp index 80e3ef7b..08edc8ae 100644 --- a/opm/core/props/BlackoilPropertiesInterface.hpp +++ b/opm/core/props/BlackoilPropertiesInterface.hpp @@ -47,6 +47,10 @@ namespace Opm /// \return N, the number of cells. virtual int numCells() const = 0; + /// Return an array containing the PVT table index for each + /// grid cell + virtual const int* cellPvtRegionIndex() const = 0; + /// \return Array of N porosity values. virtual const double* porosity() const = 0; @@ -102,14 +106,16 @@ namespace Opm /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices /// are assumed to be in Fortran order, and are typically the result /// of a call to the method matrix(). + /// \param[in] cells The index of the grid cell of each data point. /// \param[out] rho Array of nP density values, array must be valid before calling. virtual void density(const int n, const double* A, + const int* cells, double* rho) const = 0; /// Densities of stock components at surface conditions. /// \return Array of P density values. - virtual const double* surfaceDensity() const = 0; + virtual const double* surfaceDensity(int regionIdx = 0) const = 0; /// \param[in] n Number of data points. /// \param[in] s Array of nP saturation values. diff --git a/opm/core/props/pvt/BlackoilPvtProperties.cpp b/opm/core/props/pvt/BlackoilPvtProperties.cpp index cb8121a6..d1026c00 100644 --- a/opm/core/props/pvt/BlackoilPvtProperties.cpp +++ b/opm/core/props/pvt/BlackoilPvtProperties.cpp @@ -43,9 +43,7 @@ namespace Opm } void BlackoilPvtProperties::init(Opm::DeckConstPtr deck, - int numSamples, - int numCompressedCells, - const int *compressedToCartesianCellIdx) + int numSamples) { phase_usage_ = phaseUsageFromDeck(deck); @@ -69,14 +67,6 @@ namespace Opm } } - // first, calculate the PVT table index for each compressed - // cell. This array is required to construct the PVT classes - // below. - Opm::extractPvtTableIndex(pvtTableIdx_, - deck, - numCompressedCells, - compressedToCartesianCellIdx); - // Resize the property objects container props_.resize(phase_usage_.num_phases); @@ -85,7 +75,7 @@ namespace Opm // if water is used, we require the presence of the "PVTW" // keyword for now... std::shared_ptr pvtw(new PvtConstCompr); - pvtw->initFromWater(deck->getKeyword("PVTW"), pvtTableIdx_); + pvtw->initFromWater(deck->getKeyword("PVTW")); props_[phase_usage_.phase_pos[Aqua]] = pvtw; } @@ -97,18 +87,18 @@ namespace Opm Opm::DeckKeywordConstPtr pvdoKeyword = deck->getKeyword("PVDO"); if (numSamples > 0) { auto splinePvt = std::shared_ptr(new PvtDeadSpline); - splinePvt->initFromOil(pvdoKeyword, pvtTableIdx_, numSamples); + splinePvt->initFromOil(pvdoKeyword, numSamples); props_[phase_usage_.phase_pos[Liquid]] = splinePvt; } else { auto deadPvt = std::shared_ptr(new PvtDead); - deadPvt->initFromOil(pvdoKeyword, pvtTableIdx_); + deadPvt->initFromOil(pvdoKeyword); props_[phase_usage_.phase_pos[Liquid]] = deadPvt; } } else if (deck->hasKeyword("PVTO")) { - props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"), pvtTableIdx_)); + props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"))); } else if (deck->hasKeyword("PVCDO")) { std::shared_ptr pvcdo(new PvtConstCompr); - pvcdo->initFromOil(deck->getKeyword("PVCDO"), pvtTableIdx_); + pvcdo->initFromOil(deck->getKeyword("PVCDO")); props_[phase_usage_.phase_pos[Liquid]] = pvcdo; } else { @@ -123,17 +113,17 @@ namespace Opm if (numSamples > 0) { std::shared_ptr splinePvt(new PvtDeadSpline); - splinePvt->initFromGas(pvdgKeyword, pvtTableIdx_, numSamples); + splinePvt->initFromGas(pvdgKeyword, numSamples); props_[phase_usage_.phase_pos[Vapour]] = splinePvt; } else { std::shared_ptr deadPvt(new PvtDead); - deadPvt->initFromGas(pvdgKeyword, pvtTableIdx_); + deadPvt->initFromGas(pvdgKeyword); props_[phase_usage_.phase_pos[Vapour]] = deadPvt; } } else if (deck->hasKeyword("PVTG")) { - props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG"), pvtTableIdx_)); + props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG"))); } else { OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n"); } @@ -168,13 +158,14 @@ namespace Opm void BlackoilPvtProperties::mu(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_mu) const { data1_.resize(n); for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { - props_[phase]->mu(n, p, z, &data1_[0]); + props_[phase]->mu(n, pvtTableIdx, p, z, &data1_[0]); // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_mu[phase_usage_.num_phases*i + phase] = data1_[i]; @@ -183,13 +174,14 @@ namespace Opm } void BlackoilPvtProperties::B(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_B) const { data1_.resize(n); for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { - props_[phase]->B(n, p, z, &data1_[0]); + props_[phase]->B(n, pvtTableIdx, p, z, &data1_[0]); // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[phase_usage_.num_phases*i + phase] = data1_[i]; @@ -198,6 +190,7 @@ namespace Opm } void BlackoilPvtProperties::dBdp(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_B, @@ -206,7 +199,7 @@ namespace Opm data1_.resize(n); data2_.resize(n); for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { - props_[phase]->dBdp(n, p, z, &data1_[0], &data2_[0]); + props_[phase]->dBdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]); // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[phase_usage_.num_phases*i + phase] = data1_[i]; @@ -217,13 +210,14 @@ namespace Opm void BlackoilPvtProperties::R(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_R) const { data1_.resize(n); for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { - props_[phase]->R(n, p, z, &data1_[0]); + props_[phase]->R(n, pvtTableIdx, p, z, &data1_[0]); // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_R[phase_usage_.num_phases*i + phase] = data1_[i]; @@ -232,6 +226,7 @@ namespace Opm } void BlackoilPvtProperties::dRdp(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_R, @@ -240,7 +235,7 @@ namespace Opm data1_.resize(n); data2_.resize(n); for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { - props_[phase]->dRdp(n, p, z, &data1_[0], &data2_[0]); + props_[phase]->dRdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]); // #pragma omp parallel for for (int i = 0; i < n; ++i) { output_R[phase_usage_.num_phases*i + phase] = data1_[i]; diff --git a/opm/core/props/pvt/BlackoilPvtProperties.hpp b/opm/core/props/pvt/BlackoilPvtProperties.hpp index 82921c50..eea4e953 100644 --- a/opm/core/props/pvt/BlackoilPvtProperties.hpp +++ b/opm/core/props/pvt/BlackoilPvtProperties.hpp @@ -33,10 +33,12 @@ namespace Opm { /// Class collecting the pvt properties for all active phases. - /// For all the methods, the following apply: p and z - /// are expected to be of size n and n*num_phases, respectively. - /// Output arrays shall be of size n*num_phases, and must be valid - /// before calling the method. + /// For all the methods, the following apply: + /// - p and z are expected to be of size n and n*num_phases, respectively. + /// - pvtTableIdx specifies the PVT table to be used for each data + /// point and is thus expected to be an array of size n + /// - Output arrays shall be of size n*num_phases, and must be valid + /// before calling the method. /// NOTE: The difference between this interface and the one defined /// by PvtInterface is that this collects all phases' properties, /// and therefore the output arrays are of size n*num_phases as opposed @@ -49,16 +51,9 @@ namespace Opm /// Initialize from deck. /// - /// This method needs the mapping for compressed to cartesian - /// cell indices because the methods which do the actual work - /// are specified for compressed cells, but the eclipse input - /// data is specified on cartesian cell indices... - /// /// \param deck An input deck from the opm-parser module. void init(Opm::DeckConstPtr deck, - int samples, - int numCompressedCells, - const int *compressedToCartesianCellIdx); + int samples); /// \return Object describing the active phases. PhaseUsage phaseUsage() const; @@ -82,18 +77,21 @@ namespace Opm /// Viscosity as a function of p and z. void mu(const int n, + const int *pvtTableIdx, const double* p, const double* z, double* output_mu) const; /// Formation volume factor as a function of p and z. void B(const int n, + const int *pvtTableIdx, const double* p, const double* z, double* output_B) const; /// Formation volume factor and p-derivative as functions of p and z. void dBdp(const int n, + const int *pvtTableIdx, const double* p, const double* z, double* output_B, @@ -101,12 +99,14 @@ namespace Opm /// Solution factor as a function of p and z. void R(const int n, + const int *pvtTableIdx, const double* p, const double* z, double* output_R) const; /// Solution factor and p-derivative as functions of p and z. void dRdp(const int n, + const int *pvtTableIdx, const double* p, const double* z, double* output_R, @@ -119,19 +119,11 @@ namespace Opm PhaseUsage phase_usage_; - // The PVT properties. One object per active fluid phase. + // The PVT properties. We need to store one object per PVT + // region per active fluid phase. std::vector > props_; - - // The index of the PVT table which ought to be used for each - // cell. Eclipse does not seem to allow specifying fluid-phase - // specific table indices, so for the sake of simplicity, we - // don't do this either. (if it turns out that Eclipes does in - // fact support it or if it by some miracle gains this feature - // in the future, this attribute needs to become a vector of - // vectors of ints.) - std::vector pvtTableIdx_; - std::vector > densities_; + mutable std::vector data1_; mutable std::vector data2_; }; diff --git a/opm/core/props/pvt/PvtConstCompr.hpp b/opm/core/props/pvt/PvtConstCompr.hpp index e7606a06..faadd6fb 100644 --- a/opm/core/props/pvt/PvtConstCompr.hpp +++ b/opm/core/props/pvt/PvtConstCompr.hpp @@ -34,22 +34,24 @@ namespace Opm { /// Class for constant compressible phases (PVTW or PVCDO). - /// The PVT properties can either be given as a function of pressure (p) and surface volume (z) - /// or pressure (p) and gas resolution factor (r). - /// For all the virtual methods, the following apply: p, r and z - /// are expected to be of size n, size n and n*num_phases, respectively. - /// Output arrays shall be of size n, and must be valid before - /// calling the method. + /// The PVT properties can either be given as a function of + /// pressure (p) and surface volume (z) or pressure (p) and gas + /// resolution factor (r). Also, since this class supports + /// multiple PVT regions, the concrete table to be used for each + /// data point needs to be specified via the pvtTableIdx argument + /// of the respective method. For all the virtual methods, the + /// following apply: pvtTableIdx, p, r and z are expected to be of + /// size n, size n, size n and n*num_phases, respectively. Output + /// arrays shall be of size n, and must be valid before calling + /// the method. class PvtConstCompr : public PvtInterface { public: PvtConstCompr() {} - void initFromWater(Opm::DeckKeywordConstPtr pvtwKeyword, const std::vector &pvtTableIdx) + void initFromWater(Opm::DeckKeywordConstPtr pvtwKeyword) { - pvtTableIdx_ = pvtTableIdx; - int numRegions = pvtwKeyword->size(); ref_press_.resize(numRegions); @@ -69,10 +71,8 @@ namespace Opm } } - void initFromOil(Opm::DeckKeywordConstPtr pvcdoKeyword, const std::vector &pvtTableIdx) + void initFromOil(Opm::DeckKeywordConstPtr pvcdoKeyword) { - pvtTableIdx_ = pvtTableIdx; - int numRegions = pvcdoKeyword->size(); ref_press_.resize(numRegions); @@ -92,7 +92,11 @@ namespace Opm } } - PvtConstCompr(double visc) + /*! + * \brief Create a PVT object with a given viscosity that + * assumes all fluid phases to be incompressible. + */ + explicit PvtConstCompr(double visc) : ref_press_(1, 0.0), ref_B_(1, 1.0), comp_(1, 0.0), @@ -106,6 +110,7 @@ namespace Opm } virtual void mu(const int n, + const int* pvtRegionIdx, const double* p, const double* /*z*/, double* output_mu) const @@ -113,13 +118,14 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(i); + int tableIdx = getTableIndex_(pvtRegionIdx, i); double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); output_mu[i] = viscosity_[tableIdx]/(1.0 + x + 0.5*x*x); } } virtual void mu(const int n, + const int* pvtRegionIdx, const double* p, const double* /*r*/, double* output_mu, @@ -129,7 +135,7 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(i); + int tableIdx = getTableIndex_(pvtRegionIdx, i); double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); double d = (1.0 + x + 0.5*x*x); output_mu[i] = viscosity_[tableIdx]/d; @@ -139,6 +145,7 @@ namespace Opm } virtual void mu(const int n, + const int* pvtRegionIdx, const double* p, const double* /*r*/, const PhasePresence* /*cond*/, @@ -149,7 +156,7 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(i); + int tableIdx = getTableIndex_(pvtRegionIdx, i); double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); double d = (1.0 + x + 0.5*x*x); output_mu[i] = viscosity_[tableIdx]/d; @@ -159,6 +166,7 @@ namespace Opm } virtual void B(const int n, + const int* pvtRegionIdx, const double* p, const double* /*z*/, double* output_B) const @@ -166,13 +174,14 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(i); + int tableIdx = getTableIndex_(pvtRegionIdx, i); double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); output_B[i] = ref_B_[tableIdx]/(1.0 + x + 0.5*x*x); } } virtual void dBdp(const int n, + const int* pvtRegionIdx, const double* p, const double* /*z*/, double* output_B, @@ -180,7 +189,7 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int tableIdx = getTableIndex_(i); + int tableIdx = getTableIndex_(pvtRegionIdx, i); double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); double d = (1.0 + x + 0.5*x*x); output_B[i] = ref_B_[tableIdx]/d; @@ -189,6 +198,7 @@ namespace Opm } virtual void b(const int n, + const int* pvtRegionIdx, const double* p, const double* /*r*/, double* output_b, @@ -198,7 +208,7 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(i); + int tableIdx = getTableIndex_(pvtRegionIdx, i); double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); double d = (1.0 + x + 0.5*x*x); @@ -211,6 +221,7 @@ namespace Opm } virtual void b(const int n, + const int* pvtRegionIdx, const double* p, const double* /*r*/, const PhasePresence* /*cond*/, @@ -221,7 +232,7 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(i); + int tableIdx = getTableIndex_(pvtRegionIdx, i); double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); double d = (1.0 + x + 0.5*x*x); @@ -233,6 +244,7 @@ namespace Opm } virtual void rsSat(const int n, + const int* /*pvtRegionIdx*/, const double* /*p*/, double* output_rsSat, double* output_drsSatdp) const @@ -242,6 +254,7 @@ namespace Opm } virtual void rvSat(const int n, + const int* /*pvtRegionIdx*/, const double* /*p*/, double* output_rvSat, double* output_drvSatdp) const @@ -251,6 +264,7 @@ namespace Opm } virtual void R(const int n, + const int* /*pvtRegionIdx*/, const double* /*p*/, const double* /*z*/, double* output_R) const @@ -259,6 +273,7 @@ namespace Opm } virtual void dRdp(const int n, + const int* /*pvtRegionIdx*/, const double* /*p*/, const double* /*z*/, double* output_R, @@ -269,14 +284,15 @@ namespace Opm } private: - int getTableIndex_(int cellIdx) const + int getTableIndex_(const int* pvtTableIdx, int cellIdx) const { - if (pvtTableIdx_.empty()) + if (!pvtTableIdx) return 0; - return pvtTableIdx_[cellIdx]; + return pvtTableIdx[cellIdx]; } - std::vector pvtTableIdx_; + // The PVT properties. We need to store one value per PVT + // region. std::vector ref_press_; std::vector ref_B_; std::vector comp_; diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp index ecd4e885..09dc32d3 100644 --- a/opm/core/props/pvt/PvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -34,11 +34,8 @@ namespace Opm // Member functions //------------------------------------------------------------------------- /// Constructor - void PvtDead::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword, - const std::vector &pvtTableIdx) + void PvtDead::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword) { - pvtTableIdx_ = pvtTableIdx; - int numRegions = Opm::PvdoTable::numTables(pvdoKeyword); // resize the attributes of the object @@ -63,11 +60,8 @@ namespace Opm } - void PvtDead::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword, - const std::vector &pvtTableIdx) + void PvtDead::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword) { - pvtTableIdx_ = pvtTableIdx; - int numRegions = Opm::PvdgTable::numTables(pvdgKeyword); // resize the attributes of the object @@ -100,18 +94,20 @@ namespace Opm void PvtDead::mu(const int n, + const int* pvtTableIdx, const double* p, const double* /*z*/, double* output_mu) const { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_mu[i] = viscosity_[regionIdx](p[i]); } } void PvtDead::mu(const int n, + const int* pvtTableIdx, const double* p, const double* /*r*/, double* output_mu, @@ -120,7 +116,7 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_mu[i] = viscosity_[regionIdx](p[i]); output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); } @@ -129,6 +125,7 @@ namespace Opm } void PvtDead::mu(const int n, + const int* pvtTableIdx, const double* p, const double* /*r*/, const PhasePresence* /*cond*/, @@ -138,7 +135,7 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_mu[i] = viscosity_[regionIdx](p[i]); output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); } @@ -147,6 +144,7 @@ namespace Opm } void PvtDead::B(const int n, + const int* pvtTableIdx, const double* p, const double* /*z*/, double* output_B) const @@ -154,27 +152,29 @@ namespace Opm // #pragma omp parallel for // B = 1/b for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_B[i] = 1.0/b_[regionIdx](p[i]); } } void PvtDead::dBdp(const int n, + const int* pvtTableIdx, const double* p, const double* /*z*/, double* output_B, double* output_dBdp) const { - B(n, p, 0, output_B); + B(n, pvtTableIdx, p, 0, output_B); // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); double Bg = output_B[i]; output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]); } } void PvtDead::b(const int n, + const int* pvtTableIdx, const double* p, const double* /*r*/, double* output_b, @@ -184,7 +184,7 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_b[i] = b_[regionIdx](p[i]); output_dbdp[i] = b_[regionIdx].derivative(p[i]); @@ -195,6 +195,7 @@ namespace Opm } void PvtDead::b(const int n, + const int* pvtTableIdx, const double* p, const double* /*r*/, const PhasePresence* /*cond*/, @@ -205,7 +206,7 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_b[i] = b_[regionIdx](p[i]); output_dbdp[i] = b_[regionIdx].derivative(p[i]); @@ -216,6 +217,7 @@ namespace Opm } void PvtDead::rsSat(const int n, + const int* pvtTableIdx, const double* /*p*/, double* output_rsSat, double* output_drsSatdp) const @@ -225,6 +227,7 @@ namespace Opm } void PvtDead::rvSat(const int n, + const int* pvtTableIdx, const double* /*p*/, double* output_rvSat, double* output_drvSatdp) const @@ -234,6 +237,7 @@ namespace Opm } void PvtDead::R(const int n, + const int* pvtTableIdx, const double* /*p*/, const double* /*z*/, double* output_R) const @@ -242,6 +246,7 @@ namespace Opm } void PvtDead::dRdp(const int n, + const int* pvtTableIdx, const double* /*p*/, const double* /*z*/, double* output_R, diff --git a/opm/core/props/pvt/PvtDead.hpp b/opm/core/props/pvt/PvtDead.hpp index c7d26e6b..3d01eb2a 100644 --- a/opm/core/props/pvt/PvtDead.hpp +++ b/opm/core/props/pvt/PvtDead.hpp @@ -33,25 +33,28 @@ namespace Opm { /// Class for immiscible dead oil and dry gas. - /// The PVT properties can either be given as a function of pressure (p) and surface volume (z) - /// or pressure (p) and gas resolution factor (r). - /// For all the virtual methods, the following apply: p, r and z - /// are expected to be of size n, size n and n*num_phases, respectively. - /// Output arrays shall be of size n, and must be valid before - /// calling the method. + /// The PVT properties can either be given as a function of + /// pressure (p) and surface volume (z) or pressure (p) and gas + /// resolution factor (r). Also, since this class supports + /// multiple PVT regions, the concrete table to be used for each + /// data point needs to be specified via the pvtTableIdx argument + /// of the respective method. For all the virtual methods, the + /// following apply: pvtTableIdx, p, r and z are expected to be of + /// size n, size n, size n and n*num_phases, respectively. Output + /// arrays shall be of size n, and must be valid before calling + /// the method. class PvtDead : public PvtInterface { public: PvtDead() {}; - void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword, - const std::vector &pvtTableIdx); - void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword, - const std::vector &pvtTableIdx); + void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword); + void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword); virtual ~PvtDead(); /// Viscosity as a function of p and z. virtual void mu(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_mu) const; @@ -59,6 +62,7 @@ namespace Opm /// Viscosity and its derivatives as a function of p and r. /// The fluid is considered saturated if r >= rsSat(p). virtual void mu(const int n, + const int* pvtTableIdx, const double* p, const double* r, double* output_mu, @@ -68,6 +72,7 @@ namespace Opm /// Viscosity as a function of p and r. /// State condition determined by 'cond'. virtual void mu(const int n, + const int* pvtTableIdx, const double* p, const double* r, const PhasePresence* cond, @@ -77,12 +82,14 @@ namespace Opm /// Formation volume factor as a function of p and z. virtual void B(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_B) const; /// Formation volume factor and p-derivative as functions of p and z. virtual void dBdp(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_B, @@ -91,6 +98,7 @@ namespace Opm /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The fluid is considered saturated if r >= rsSat(p). virtual void b(const int n, + const int* pvtTableIdx, const double* p, const double* r, double* output_b, @@ -100,6 +108,7 @@ namespace Opm /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// State condition determined by 'cond'. virtual void b(const int n, + const int* pvtTableIdx, const double* p, const double* r, const PhasePresence* cond, @@ -110,31 +119,42 @@ namespace Opm /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p. virtual void rsSat(const int n, + const int* pvtTableIdx, const double* p, double* output_rsSat, double* output_drsSatdp) const; /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p. virtual void rvSat(const int n, + const int* pvtTableIdx, const double* p, double* output_rvSat, double* output_drvSatdp) const; /// Solution factor as a function of p and z. virtual void R(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_R) const; /// Solution factor and p-derivative as functions of p and z. virtual void dRdp(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_R, double* output_dRdp) const; private: - // PVT properties of dry gas or dead oil - std::vector pvtTableIdx_; + int getTableIndex_(const int* pvtTableIdx, int cellIdx) const + { + if (!pvtTableIdx) + return 0; + return pvtTableIdx[cellIdx]; + } + + // PVT properties of dry gas or dead oil. We need to store one + // table per PVT region. std::vector > b_; std::vector > viscosity_; }; diff --git a/opm/core/props/pvt/PvtDeadSpline.cpp b/opm/core/props/pvt/PvtDeadSpline.cpp index 1b40fae4..69f749b7 100644 --- a/opm/core/props/pvt/PvtDeadSpline.cpp +++ b/opm/core/props/pvt/PvtDeadSpline.cpp @@ -40,11 +40,8 @@ namespace Opm {} void PvtDeadSpline::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword, - const std::vector &pvtTableIdx, int numSamples) { - pvtTableIdx_ = pvtTableIdx; - int numRegions = Opm::PvdoTable::numTables(pvdoKeyword); // resize the attributes of the object @@ -72,11 +69,8 @@ namespace Opm } void PvtDeadSpline::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword, - const std::vector &pvtTableIdx, int numSamples) { - pvtTableIdx_ = pvtTableIdx; - int numRegions = Opm::PvdgTable::numTables(pvdgKeyword); // resize the attributes of the object @@ -111,18 +105,20 @@ namespace Opm void PvtDeadSpline::mu(const int n, + const int* pvtTableIdx, const double* p, const double* /*z*/, double* output_mu) const { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_mu[i] = viscosity_[regionIdx](p[i]); } } void PvtDeadSpline::mu(const int n, + const int* pvtTableIdx, const double* p, const double* /*r*/, double* output_mu, @@ -131,7 +127,7 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_mu[i] = viscosity_[regionIdx](p[i]); output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); } @@ -139,6 +135,7 @@ namespace Opm } void PvtDeadSpline::mu(const int n, + const int* pvtTableIdx, const double* p, const double* /*r*/, const PhasePresence* /*cond*/, @@ -149,7 +146,7 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_mu[i] = viscosity_[regionIdx](p[i]); output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); } @@ -157,33 +154,36 @@ namespace Opm } void PvtDeadSpline::B(const int n, + const int* pvtTableIdx, const double* p, const double* /*z*/, double* output_B) const { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_B[i] = 1.0/b_[regionIdx](p[i]); } } void PvtDeadSpline::dBdp(const int n, + const int* pvtTableIdx, const double* p, const double* /*z*/, double* output_B, double* output_dBdp) const { - B(n, p, 0, output_B); + B(n, pvtTableIdx, p, 0, output_B); // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); double Bg = output_B[i]; output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]); } } void PvtDeadSpline::b(const int n, + const int* pvtTableIdx, const double* p, const double* /*r*/, double* output_b, @@ -192,7 +192,7 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_b[i] = b_[regionIdx](p[i]); output_dbdp[i] = b_[regionIdx].derivative(p[i]); } @@ -200,6 +200,7 @@ namespace Opm } void PvtDeadSpline::b(const int n, + const int* pvtTableIdx, const double* p, const double* /*r*/, const PhasePresence* /*cond*/, @@ -209,7 +210,7 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int regionIdx = pvtTableIdx_[i]; + int regionIdx = getTableIndex_(pvtTableIdx, i); output_b[i] = b_[regionIdx](p[i]); output_dbdp[i] = b_[regionIdx].derivative(p[i]); } @@ -217,6 +218,7 @@ namespace Opm } void PvtDeadSpline::rsSat(const int n, + const int* /*pvtTableIdx*/, const double* /*p*/, double* output_rsSat, double* output_drsSatdp) const @@ -226,6 +228,7 @@ namespace Opm } void PvtDeadSpline::rvSat(const int n, + const int* /*pvtTableIdx*/, const double* /*p*/, double* output_rvSat, double* output_drvSatdp) const @@ -235,6 +238,7 @@ namespace Opm } void PvtDeadSpline::R(const int n, + const int* /*pvtTableIdx*/, const double* /*p*/, const double* /*z*/, double* output_R) const @@ -243,6 +247,7 @@ namespace Opm } void PvtDeadSpline::dRdp(const int n, + const int* /*pvtTableIdx*/, const double* /*p*/, const double* /*z*/, double* output_R, diff --git a/opm/core/props/pvt/PvtDeadSpline.hpp b/opm/core/props/pvt/PvtDeadSpline.hpp index b47a4986..b7c8beba 100644 --- a/opm/core/props/pvt/PvtDeadSpline.hpp +++ b/opm/core/props/pvt/PvtDeadSpline.hpp @@ -44,16 +44,15 @@ namespace Opm PvtDeadSpline(); void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword, - const std::vector &pvtTableIdx, int numSamples); void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword, - const std::vector &pvtTableIdx, int numSamples); virtual ~PvtDeadSpline(); /// Viscosity as a function of p and z. virtual void mu(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_mu) const; @@ -61,6 +60,7 @@ namespace Opm /// Viscosity and its derivatives as a function of p and r. /// The fluid is considered saturated if r >= rsSat(p). virtual void mu(const int n, + const int* pvtTableIdx, const double* p, const double* r, double* output_mu, @@ -70,6 +70,7 @@ namespace Opm /// Viscosity as a function of p and r. /// State condition determined by 'cond'. virtual void mu(const int n, + const int* pvtTableIdx, const double* p, const double* r, const PhasePresence* cond, @@ -79,12 +80,14 @@ namespace Opm /// Formation volume factor as a function of p and z. virtual void B(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_B) const; /// Formation volume factor and p-derivative as functions of p and z. virtual void dBdp(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_B, @@ -93,6 +96,7 @@ namespace Opm /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The fluid is considered saturated if r >= rsSat(p). virtual void b(const int n, + const int* pvtTableIdx, const double* p, const double* r, double* output_b, @@ -102,6 +106,7 @@ namespace Opm /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// State condition determined by 'cond'. virtual void b(const int n, + const int* pvtTableIdx, const double* p, const double* r, const PhasePresence* cond, @@ -111,31 +116,42 @@ namespace Opm /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p. virtual void rsSat(const int n, + const int* pvtTableIdx, const double* p, double* output_rsSat, double* output_drsSatdp) const; /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p. virtual void rvSat(const int n, + const int* pvtTableIdx, const double* p, double* output_rvSat, double* output_drvSatdp) const; /// Solution factor as a function of p and z. virtual void R(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_R) const; /// Solution factor and p-derivative as functions of p and z. virtual void dRdp(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_R, double* output_dRdp) const; private: - // PVT properties of dry gas or dead oil - std::vector pvtTableIdx_; + int getTableIndex_(const int* pvtTableIdx, int cellIdx) const + { + if (!pvtTableIdx) + return 0; + return pvtTableIdx[cellIdx]; + } + + // PVT properties of dry gas or dead oil. We need to store one + // table per PVT region. std::vector > b_; std::vector > viscosity_; }; diff --git a/opm/core/props/pvt/PvtInterface.hpp b/opm/core/props/pvt/PvtInterface.hpp index 8be7437a..85b5b74f 100644 --- a/opm/core/props/pvt/PvtInterface.hpp +++ b/opm/core/props/pvt/PvtInterface.hpp @@ -44,13 +44,19 @@ namespace Opm /// The PVT properties can either be given as a function of pressure (p) and surface volume (z) /// or pressure (p) and gas resolution factor (r). - /// For all the virtual methods, the following apply: p, r and z - /// are expected to be of size n, size n and n*num_phases, respectively. - /// Output arrays shall be of size n, and must be valid before - /// calling the method. + /// For all the virtual methods, the following apply: + /// - pvtRegionIdx is an array of size n and represents the + /// index of the PVT table which should be used to calculate + /// the output. NULL can also be passed and is interpreted + /// such that the first table should be used for the output + /// - p, r and z are expected to be of size n, size n and + /// n*num_phases, respectively. + /// - Output arrays shall be of size n, and must be valid before + /// calling the method. /// Viscosity as a function of p and z. virtual void mu(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_mu) const = 0; @@ -58,6 +64,7 @@ namespace Opm /// Viscosity as a function of p and r. /// The fluid is considered saturated if r >= rsSat(p). virtual void mu(const int n, + const int* pvtRegionIdx, const double* p, const double* r, double* output_mu, @@ -67,6 +74,7 @@ namespace Opm /// Viscosity as a function of p and r. /// State condition determined by 'cond'. virtual void mu(const int n, + const int* pvtRegionIdx, const double* p, const double* r, const PhasePresence* cond, @@ -76,12 +84,14 @@ namespace Opm /// Formation volume factor as a function of p and z. virtual void B(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_B) const = 0; /// Formation volume factor and p-derivative as functions of p and z. virtual void dBdp(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_B, @@ -90,6 +100,7 @@ namespace Opm /// The inverse of the volume factor b = 1 / B as a function of p and r. /// The fluid is considered saturated if r >= rsSat(p). virtual void b(const int n, + const int* pvtRegionIdx, const double* p, const double* r, double* output_b, @@ -99,6 +110,7 @@ namespace Opm /// The inverse of the volume factor b = 1 / B as a function of p and r. /// State condition determined by 'cond'. virtual void b(const int n, + const int* pvtRegionIdx, const double* p, const double* r, const PhasePresence* cond, @@ -108,12 +120,14 @@ namespace Opm /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p. virtual void rsSat(const int n, + const int* pvtRegionIdx, const double* p, double* output_rsSat, double* output_drsSatdp) const = 0; /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p. virtual void rvSat(const int n, + const int* pvtRegionIdx, const double* p, double* output_rvSat, double* output_drvSatdp) const = 0; @@ -121,12 +135,14 @@ namespace Opm /// Solution factor as a function of p and z. virtual void R(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_R) const = 0; /// Solution factor and p-derivative as functions of p and z. virtual void dRdp(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_R, @@ -140,18 +156,21 @@ namespace Opm /*! * \brief Helper function to create an array containing the (C-Style) - * PVT table index for each compressed cell. + * PVT table index for each compressed cell from an Eclipse deck. * - * The main point of this function is to avoid code duplication - * because the Eclipse deck only contains Fortran-style PVT table - * indices which start at 1 instead of 0 and -- more relevantly -- it - * uses logically cartesian cell indices to specify the table index of - * a cell. + * This function assumes that the degrees of freedom where PVT + * properties need to be calculated are grid cells. The main point + * of this function is to avoid code duplication because the + * Eclipse deck only contains Fortran-style PVT table indices + * which start at 1 instead of 0 and -- more relevantly -- it uses + * logically cartesian cell indices to specify the table index of + * a cell while the classes which use the PvtInterface + * implementations usually use compressed cells. */ - void extractPvtTableIndex(std::vector &pvtTableIdx, + void extractPvtTableIndex(std::vector& pvtTableIdx, Opm::DeckConstPtr deck, size_t numCompressed, - const int *compressedToCartesianIdx); + const int* compressedToCartesianIdx); } // namespace Opm diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp index 6e1109a6..7c59cc3a 100644 --- a/opm/core/props/pvt/PvtLiveGas.cpp +++ b/opm/core/props/pvt/PvtLiveGas.cpp @@ -46,11 +46,8 @@ namespace Opm //------------------------------------------------------------------------ // Member functions //------------------------------------------------------------------------- - PvtLiveGas::PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword, - const std::vector &pvtTableIdx) + PvtLiveGas::PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword) { - pvtTableIdx_ = pvtTableIdx; - int numTables = Opm::PvtgTable::numTables(pvtgKeyword); saturated_gas_table_.resize(numTables); undersat_gas_tables_.resize(numTables); @@ -93,18 +90,20 @@ namespace Opm void PvtLiveGas::mu(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_mu) const { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - output_mu[i] = miscible_gas(p[i], z + num_phases_*i, pvtTableIdx_[i], 2, false); + output_mu[i] = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 2, false); } } /// Viscosity and its derivatives as a function of p and r. void PvtLiveGas::mu(const int /*n*/, + const int* /*pvtRegionIdx*/, const double* /*p*/, const double* /*r*/, double* /*output_mu*/, @@ -116,6 +115,7 @@ namespace Opm /// Viscosity and its derivatives as a function of p and r. void PvtLiveGas::mu(const int n, + const int* pvtRegionIdx, const double* p, const double* r, const PhasePresence* cond, @@ -125,9 +125,10 @@ namespace Opm { for (int i = 0; i < n; ++i) { const PhasePresence& cnd = cond[i]; - output_mu[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 2, 0); - output_dmudp[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 2, 1); - output_dmudr[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 2, 2); + 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); } } @@ -135,13 +136,14 @@ namespace Opm /// Formation volume factor as a function of p and z. void PvtLiveGas::B(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_B) const { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - output_B[i] = evalB(p[i], z + num_phases_*i, pvtTableIdx_[i]); + output_B[i] = evalB(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i)); } } @@ -149,6 +151,7 @@ namespace Opm /// Formation volume factor and p-derivative as functions of p and z. void PvtLiveGas::dBdp(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_B, @@ -156,12 +159,13 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - evalBDeriv(p[i], z + num_phases_*i, pvtTableIdx_[i], output_B[i], output_dBdp[i]); + evalBDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_B[i], output_dBdp[i]); } } /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. void PvtLiveGas::b(const int /*n*/, + const int* /*pvtRegionIdx*/, const double* /*p*/, const double* /*r*/, double* /*output_b*/, @@ -174,6 +178,7 @@ namespace Opm /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. void PvtLiveGas::b(const int n, + const int* pvtRegionIdx, const double* p, const double* r, const PhasePresence* cond, @@ -186,21 +191,23 @@ namespace Opm for (int i = 0; i < n; ++i) { const PhasePresence& cnd = cond[i]; - output_b[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 1, 0); - output_dbdp[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 1, 1); - output_dbdr[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 1, 2); + int tableIdx = getTableIndex_(pvtRegionIdx, i); + output_b[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 0); + output_dbdp[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 1); + output_dbdr[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 2); } } /// Gas resolution and its derivatives at bublepoint as a function of p. void PvtLiveGas::rvSat(const int n, + const int* pvtRegionIdx, const double* p, double* output_rvSat, double* output_drvSatdp) const { for (int i = 0; i < n; ++i) { - int pvtTableIdx = pvtTableIdx_[i]; + int pvtTableIdx = getTableIndex_(pvtRegionIdx, i); output_rvSat[i] = linearInterpolation(saturated_gas_table_[pvtTableIdx][0], saturated_gas_table_[pvtTableIdx][3],p[i]); output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0], @@ -210,6 +217,7 @@ namespace Opm } void PvtLiveGas::rsSat(const int n, + const int* pvtRegionIdx, const double* /*p*/, double* output_rsSat, double* output_drsSatdp) const @@ -220,13 +228,14 @@ namespace Opm /// Solution factor as a function of p and z. void PvtLiveGas::R(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_R) const { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - output_R[i] = evalR(p[i], z + num_phases_*i, pvtTableIdx_[i]); + output_R[i] = evalR(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i)); } } @@ -234,6 +243,7 @@ namespace Opm /// Solution factor and p-derivative as functions of p and z. void PvtLiveGas::dRdp(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_R, @@ -241,7 +251,7 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - evalRDeriv(p[i], z + num_phases_*i, pvtTableIdx_[i], output_R[i], output_dRdp[i]); + evalRDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_R[i], output_dRdp[i]); } } diff --git a/opm/core/props/pvt/PvtLiveGas.hpp b/opm/core/props/pvt/PvtLiveGas.hpp index 7cdc361e..86ed258e 100644 --- a/opm/core/props/pvt/PvtLiveGas.hpp +++ b/opm/core/props/pvt/PvtLiveGas.hpp @@ -38,11 +38,12 @@ namespace Opm class PvtLiveGas : public PvtInterface { public: - PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword, const std::vector &pvtTableIdx); + PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword); virtual ~PvtLiveGas(); /// Viscosity as a function of p and z. virtual void mu(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_mu) const; @@ -50,6 +51,7 @@ namespace Opm /// Viscosity and its derivatives as a function of p and r. /// The fluid is considered saturated if r >= rsSat(p). virtual void mu(const int n, + const int* pvtRegionIdx, const double* p, const double* r, double* output_mu, @@ -59,6 +61,7 @@ namespace Opm /// Viscosity as a function of p and r. /// State condition determined by 'cond'. virtual void mu(const int n, + const int* pvtRegionIdx, const double* p, const double* r, const PhasePresence* cond, @@ -68,12 +71,14 @@ namespace Opm /// Formation volume factor as a function of p and z. virtual void B(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_B) const; /// Formation volume factor and p-derivative as functions of p and z. virtual void dBdp(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_B, @@ -82,6 +87,7 @@ namespace Opm /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The fluid is considered saturated if r >= rsSat(p). virtual void b(const int n, + const int* pvtRegionIdx, const double* p, const double* r, double* output_b, @@ -91,6 +97,7 @@ namespace Opm /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// State condition determined by 'cond'. virtual void b(const int n, + const int* pvtRegionIdx, const double* p, const double* r, const PhasePresence* cond, @@ -102,30 +109,41 @@ namespace Opm /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p. virtual void rsSat(const int n, + const int* pvtRegionIdx, const double* p, double* output_rsSat, double* output_drsSatdp) const; /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p. virtual void rvSat(const int n, + const int* pvtRegionIdx, const double* p, double* output_rvSat, double* output_drvSatdp) const; /// Solution factor as a function of p and z. virtual void R(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_R) const; /// Solution factor and p-derivative as functions of p and z. virtual void dRdp(const int n, + const int* pvtRegionIdx, const double* p, const double* z, double* output_R, double* output_dRdp) const; protected: + int getTableIndex_(const int* pvtTableIdx, int cellIdx) const + { + if (!pvtTableIdx) + return 0; + return pvtTableIdx[cellIdx]; + } + double evalB(double press, const double* surfvol, int pvtTableIdx) const; void evalBDeriv(double press, const double* surfvol, int pvtTableIdx, double& B, double& dBdp) const; double evalR(double press, const double* surfvol, int pvtTableIdx) const; @@ -143,8 +161,8 @@ namespace Opm const int pvtTableIdx, const int item, const int deriv = 0) const; - // PVT properties of wet gas (with vaporised oil) - std::vector pvtTableIdx_; + // PVT properties of wet gas (with vaporised oil). We need to + // store one table per PVT region. std::vector< std::vector > > saturated_gas_table_; std::vector< std::vector > > > undersat_gas_tables_; }; diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp index 2c6691f8..578e6369 100644 --- a/opm/core/props/pvt/PvtLiveOil.cpp +++ b/opm/core/props/pvt/PvtLiveOil.cpp @@ -34,11 +34,8 @@ namespace Opm //------------------------------------------------------------------------ // Member functions //------------------------------------------------------------------------- - PvtLiveOil::PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword, - const std::vector &pvtTableIdx) + PvtLiveOil::PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword) { - pvtTableIdx_ = pvtTableIdx; - int numTables = Opm::PvtoTable::numTables(pvtoKeyword); saturated_oil_table_.resize(numTables); undersat_oil_tables_.resize(numTables); @@ -120,20 +117,22 @@ namespace Opm /// Viscosity as a function of p and z. void PvtLiveOil::mu(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_mu) const { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int pvtTableIdx = pvtTableIdx_[i]; + int tableIdx = getTableIndex_(pvtTableIdx, i); - output_mu[i] = miscible_oil(p[i], z + num_phases_*i, pvtTableIdx, 2, false); + output_mu[i] = miscible_oil(p[i], z + num_phases_*i, tableIdx, 2, false); } } /// Viscosity and its derivatives as a function of p and r. void PvtLiveOil::mu(const int n, + const int* pvtTableIdx, const double* p, const double* r, double* output_mu, @@ -142,17 +141,18 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int pvtTableIdx = pvtTableIdx_[i]; + int tableIdx = getTableIndex_(pvtTableIdx, i); - output_mu[i] = miscible_oil(p[i], r[i], pvtTableIdx, 2, 0); - output_dmudp[i] = miscible_oil(p[i], r[i], pvtTableIdx, 2, 1); - output_dmudr[i] = miscible_oil(p[i], r[i], pvtTableIdx, 2, 2); + 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); } } /// Viscosity and its derivatives as a function of p and r. void PvtLiveOil::mu(const int n, + const int* pvtTableIdx, const double* p, const double* r, const PhasePresence* cond, @@ -162,12 +162,12 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int pvtTableIdx = pvtTableIdx_[i]; + int tableIdx = getTableIndex_(pvtTableIdx, i); const PhasePresence& cnd = cond[i]; - output_mu[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 2, 0); - output_dmudp[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 2, 1); - output_dmudr[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 2, 2); + 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); } } @@ -175,15 +175,16 @@ namespace Opm /// Formation volume factor as a function of p and z. void PvtLiveOil::B(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_B) const { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int pvtTableIdx = pvtTableIdx_[i]; + int tableIdx = getTableIndex_(pvtTableIdx, i); - output_B[i] = evalB(pvtTableIdx, p[i], z + num_phases_*i); + output_B[i] = evalB(tableIdx, p[i], z + num_phases_*i); } } @@ -191,6 +192,7 @@ namespace Opm /// Formation volume factor and p-derivative as functions of p and z. void PvtLiveOil::dBdp(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_B, @@ -198,13 +200,14 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int pvtTableIdx = pvtTableIdx_[i]; + int tableIdx = getTableIndex_(pvtTableIdx, i); - evalBDeriv(pvtTableIdx, p[i], z + num_phases_*i, output_B[i], output_dBdp[i]); + evalBDeriv(tableIdx, p[i], z + num_phases_*i, output_B[i], output_dBdp[i]); } } void PvtLiveOil::b(const int n, + const int* pvtTableIdx, const double* p, const double* r, double* output_b, @@ -214,16 +217,17 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int pvtTableIdx = pvtTableIdx_[i]; + int tableIdx = getTableIndex_(pvtTableIdx, i); - output_b[i] = miscible_oil(pvtTableIdx, p[i], r[i], 1, 0); - output_dbdp[i] = miscible_oil(pvtTableIdx, p[i], r[i], 1, 1); - output_dbdr[i] = miscible_oil(pvtTableIdx, p[i], r[i], 1, 2); + output_b[i] = miscible_oil(tableIdx, p[i], r[i], 1, 0); + output_dbdp[i] = miscible_oil(tableIdx, p[i], r[i], 1, 1); + output_dbdr[i] = miscible_oil(tableIdx, p[i], r[i], 1, 2); } } void PvtLiveOil::b(const int n, + const int* pvtTableIdx, const double* p, const double* r, const PhasePresence* cond, @@ -235,32 +239,34 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { const PhasePresence& cnd = cond[i]; - int pvtTableIdx = pvtTableIdx_[i]; + int tableIdx = getTableIndex_(pvtTableIdx, i); - output_b[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 1, 0); - output_dbdp[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 1, 1); - output_dbdr[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 1, 2); + output_b[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 0); + output_dbdp[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 1); + output_dbdr[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 2); } } void PvtLiveOil::rsSat(const int n, + const int* pvtTableIdx, const double* p, double* output_rsSat, double* output_drsSatdp) const { for (int i = 0; i < n; ++i) { - int pvtTableIdx = pvtTableIdx_[i]; - output_rsSat[i] = linearInterpolation(saturated_oil_table_[pvtTableIdx][0], - saturated_oil_table_[pvtTableIdx][3],p[i]); - output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0], - saturated_oil_table_[pvtTableIdx][3],p[i]); + int tableIdx = getTableIndex_(pvtTableIdx, i); + output_rsSat[i] = linearInterpolation(saturated_oil_table_[tableIdx][0], + saturated_oil_table_[tableIdx][3],p[i]); + output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[tableIdx][0], + saturated_oil_table_[tableIdx][3],p[i]); } } void PvtLiveOil::rvSat(const int n, + const int* pvtTableIdx, const double* /*p*/, double* output_rvSat, double* output_drvSatdp) const @@ -271,14 +277,15 @@ namespace Opm /// Solution factor as a function of p and z. void PvtLiveOil::R(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_R) const { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int pvtTableIdx = pvtTableIdx_[i]; - output_R[i] = evalR(pvtTableIdx, p[i], z + num_phases_*i); + int tableIdx = getTableIndex_(pvtTableIdx, i); + output_R[i] = evalR(tableIdx, p[i], z + num_phases_*i); } } @@ -286,6 +293,7 @@ namespace Opm /// Solution factor and p-derivative as functions of p and z. void PvtLiveOil::dRdp(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_R, @@ -293,8 +301,8 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - int pvtTableIdx = pvtTableIdx_[i]; - evalRDeriv(pvtTableIdx, p[i], z + num_phases_*i, output_R[i], output_dRdp[i]); + int tableIdx = getTableIndex_(pvtTableIdx, i); + evalRDeriv(tableIdx, p[i], z + num_phases_*i, output_R[i], output_dRdp[i]); } } diff --git a/opm/core/props/pvt/PvtLiveOil.hpp b/opm/core/props/pvt/PvtLiveOil.hpp index b16708a0..b6f61644 100644 --- a/opm/core/props/pvt/PvtLiveOil.hpp +++ b/opm/core/props/pvt/PvtLiveOil.hpp @@ -39,11 +39,12 @@ namespace Opm class PvtLiveOil : public PvtInterface { public: - PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword, const std::vector &pvtTableIdx); + PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword); virtual ~PvtLiveOil(); /// Viscosity as a function of p and z. virtual void mu(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_mu) const; @@ -51,6 +52,7 @@ namespace Opm /// Viscosity and its derivatives as a function of p and r. /// The fluid is considered saturated if r >= rsSat(p). virtual void mu(const int n, + const int* pvtTableIdx, const double* p, const double* r, double* output_mu, @@ -60,6 +62,7 @@ namespace Opm /// Viscosity as a function of p and r. /// State condition determined by 'cond'. virtual void mu(const int n, + const int* pvtTableIdx, const double* p, const double* r, const PhasePresence* cond, @@ -69,12 +72,14 @@ namespace Opm /// Formation volume factor as a function of p and z. virtual void B(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_B) const; /// Formation volume factor and p-derivative as functions of p and z. virtual void dBdp(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_B, @@ -83,6 +88,7 @@ namespace Opm /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The fluid is considered saturated if r >= rsSat(p). virtual void b(const int n, + const int* pvtTableIdx, const double* p, const double* r, double* output_b, @@ -92,6 +98,7 @@ namespace Opm /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// State condition determined by 'cond'. virtual void b(const int n, + const int* pvtTableIdx, const double* p, const double* r, const PhasePresence* cond, @@ -101,30 +108,41 @@ namespace Opm /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p. virtual void rsSat(const int n, + const int* pvtTableIdx, const double* p, double* output_rsSat, double* output_drsSatdp) const; /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p. virtual void rvSat(const int n, + const int* pvtTableIdx, const double* p, double* output_rvSat, double* output_drvSatdp) const; /// Solution factor as a function of p and z. virtual void R(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_R) const; /// Solution factor and p-derivative as functions of p and z. virtual void dRdp(const int n, + const int* pvtTableIdx, const double* p, const double* z, double* output_R, double* output_dRdp) const; private: + int getTableIndex_(const int* pvtTableIdx, int cellIdx) const + { + if (!pvtTableIdx) + return 0; + return pvtTableIdx[cellIdx]; + } + double evalB(size_t pvtTableIdx, double press, const double* surfvol) const; void evalBDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& B, double& dBdp) const; double evalR(size_t pvtTableIdx, double press, const double* surfvol) const; @@ -150,8 +168,8 @@ namespace Opm const int item, const int deriv = 0) const; - // PVT properties of live oil (with dissolved gas) - std::vector pvtTableIdx_; + // PVT properties of live oil (with dissolved gas). We need to + // store one table per PVT region. std::vector > > saturated_oil_table_; std::vector > > > undersat_oil_tables_; }; diff --git a/opm/core/props/pvt/PvtPropertiesBasic.hpp b/opm/core/props/pvt/PvtPropertiesBasic.hpp index e9198ff9..46add622 100644 --- a/opm/core/props/pvt/PvtPropertiesBasic.hpp +++ b/opm/core/props/pvt/PvtPropertiesBasic.hpp @@ -95,6 +95,8 @@ namespace Opm double* output_dRdp) const; private: + // The PVT properties. We need to store one value per PVT + // region. std::vector density_; std::vector viscosity_; std::vector formation_volume_factor_; diff --git a/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp b/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp index e23839e9..27176351 100644 --- a/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp +++ b/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp @@ -33,9 +33,9 @@ namespace Opm { } - void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck ) + void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck) { - // If we need multiple regions, this class and the Pvt* classes must change. + // So far, this class only supports a single PVT region. TODO? int region_number = 0; PhaseUsage phase_usage = phaseUsageFromDeck(deck); diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index dd90e368..e97d5077 100644 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ b/opm/core/simulator/EquilibrationHelpers.hpp @@ -139,7 +139,7 @@ namespace Opm props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp); std::vector rho(np, 0.0); - props_.density(1, &A[0], &rho[0]); + props_.density(1, &A[0], &c_[0], &rho[0]); return rho; } diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index e5f8c8ec..a4ec1ee4 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -369,8 +369,6 @@ namespace Opm const UnstructuredGrid& G , const double grav) { - typedef Miscibility::NoMixing NoMix; - for (typename RMap::RegionId r = 0, nr = reg.numRegions(); r < nr; ++r) diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index 6b1885ea..70a6248d 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -184,7 +184,7 @@ namespace Opm double A[4] = { 0.0 }; props_.matrix(1, &pressure, surfvol[phase], cells, A, 0); double rho[2] = { 0.0 }; - props_.density(1, A, rho); + props_.density(1, A, cells, rho); return rho[phase]; } }; diff --git a/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp index 4370a8ef..82f8ea80 100644 --- a/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp +++ b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp @@ -418,7 +418,7 @@ namespace Opm assert(np == 2); const int dim = grid_.dimensions; density_.resize(nc*np); - props_.density(grid_.number_of_cells, &A_[0], &density_[0]); + props_.density(grid_.number_of_cells, &A_[0], /*cellIndices=*/NULL, &density_[0]); std::fill(gravflux_.begin(), gravflux_.end(), 0.0); for (int f = 0; f < nf; ++f) { const int* c = &grid_.face_cells[2*f]; diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 2454eb73..8a0ebac8 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -285,7 +285,8 @@ namespace Opm /// @brief Computes saturation from surface volume void computeSaturation(const BlackoilPropertiesInterface& props, - BlackoilState& state){ + BlackoilState& state) + { const int np = props.numPhases(); const int nc = props.numCells(); diff --git a/tests/test_blackoilfluid.cpp b/tests/test_blackoilfluid.cpp index 8dcf3417..6f5cc3ca 100644 --- a/tests/test_blackoilfluid.cpp +++ b/tests/test_blackoilfluid.cpp @@ -40,20 +40,10 @@ using namespace std; std::vector > getProps(Opm::DeckConstPtr deck, PhaseUsage phase_usage_){ Opm::GridManager grid(deck); - const UnstructuredGrid* cgrid = grid.c_grid(); enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 }; int samples = 0; - // first, calculate the PVT table index for each compressed - // cell. This array is required to construct the PVT classes - // below. - std::vector pvtTableIdx; - Opm::extractPvtTableIndex(pvtTableIdx, - deck, - cgrid->number_of_cells, - cgrid->global_cell); - std::vector > props_; // Set the properties. props_.resize(phase_usage_.num_phases); @@ -62,7 +52,7 @@ std::vector > getProps(Opm::DeckConstPtr deck, Pha if (phase_usage_.phase_used[Aqua]) { if (deck->hasKeyword("PVTW")) { std::shared_ptr pvtw(new PvtConstCompr); - pvtw->initFromWater(deck->getKeyword("PVTW"), pvtTableIdx); + pvtw->initFromWater(deck->getKeyword("PVTW")); props_[phase_usage_.phase_pos[Aqua]] = pvtw; } else { @@ -77,18 +67,18 @@ std::vector > getProps(Opm::DeckConstPtr deck, Pha Opm::DeckKeywordConstPtr pvdoKeyword(deck->getKeyword("PVDO")); if (samples > 0) { std::shared_ptr splinePvt(new PvtDeadSpline); - splinePvt->initFromOil(pvdoKeyword, pvtTableIdx, samples); + splinePvt->initFromOil(pvdoKeyword, samples); props_[phase_usage_.phase_pos[Liquid]] = splinePvt; } else { std::shared_ptr deadPvt(new PvtDead); - deadPvt->initFromOil(pvdoKeyword, pvtTableIdx); + deadPvt->initFromOil(pvdoKeyword); props_[phase_usage_.phase_pos[Liquid]] = deadPvt; } } else if (deck->hasKeyword("PVTO")) { - props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"), pvtTableIdx)); + props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"))); } else if (deck->hasKeyword("PVCDO")) { std::shared_ptr pvcdo(new PvtConstCompr); - pvcdo->initFromOil(deck->getKeyword("PVCDO"), pvtTableIdx); + pvcdo->initFromOil(deck->getKeyword("PVCDO")); props_[phase_usage_.phase_pos[Liquid]] = pvcdo; } else { @@ -101,15 +91,15 @@ std::vector > getProps(Opm::DeckConstPtr deck, Pha Opm::DeckKeywordConstPtr pvdgKeyword(deck->getKeyword("PVDG")); if (samples > 0) { std::shared_ptr splinePvt(new PvtDeadSpline); - splinePvt->initFromGas(pvdgKeyword, pvtTableIdx, samples); + splinePvt->initFromGas(pvdgKeyword, samples); props_[phase_usage_.phase_pos[Vapour]] = splinePvt; } else { std::shared_ptr deadPvt(new PvtDead); - deadPvt->initFromGas(pvdgKeyword, pvtTableIdx); + deadPvt->initFromGas(pvdgKeyword); props_[phase_usage_.phase_pos[Vapour]] = deadPvt; } } else if (deck->hasKeyword("PVTG")) { - props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG"), pvtTableIdx)); + props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG"))); } else { OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n"); } @@ -118,7 +108,7 @@ std::vector > getProps(Opm::DeckConstPtr deck, Pha return props_; } -void testmu(const double reltol, int n, int np, std::vector p, std::vector r,std::vector z, +void testmu(const double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector r,std::vector z, std::vector > props_, std::vector condition) { std::vector mu(n); @@ -132,8 +122,8 @@ void testmu(const double reltol, int n, int np, std::vector p, std::vect // 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]); + props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]); + props_[phase]->mu(n, &pvtTableIdx[0], &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]); @@ -154,7 +144,7 @@ void testmu(const double reltol, int n, int np, std::vector p, std::vect } } -void testb(const double reltol, int n, int np, std::vector p, std::vector r,std::vector z, +void testb(const double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector r,std::vector z, std::vector > props_, std::vector condition) { // test b @@ -171,8 +161,8 @@ void testb(const double reltol, int n, int np, std::vector p, std::vecto 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]); + props_[phase]->b(n, &pvtTableIdx[0], &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]); + props_[phase]->dBdp(n, &pvtTableIdx[0], &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]); @@ -197,7 +187,7 @@ void testb(const double reltol, int n, int np, std::vector p, std::vecto } } -void testrsSat(double reltol, int n, int np, std::vector p, std::vector > props_){ +void testrsSat(double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector > props_){ // test bublepoint pressure std::vector rs(n); std::vector drsdp(n); @@ -205,7 +195,7 @@ void testrsSat(double reltol, int n, int np, std::vector p, std::vector< double drsdp_diff_u; for (int phase = 0; phase < np; ++phase) { - props_[phase] ->rsSat(n, &p[0], &rs[0], &drsdp[0]); + props_[phase] ->rsSat(n, &pvtTableIdx[0], &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]); @@ -219,7 +209,7 @@ void testrsSat(double reltol, int n, int np, std::vector p, std::vector< } } -void testrvSat(double reltol, int n, int np, std::vector p, std::vector > props_){ +void testrvSat(double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector > props_){ // test rv saturated std::vector rv(n); std::vector drvdp(n); @@ -227,7 +217,7 @@ void testrvSat(double reltol, int n, int np, std::vector p, std::vector< double drvdp_diff_u; for (int phase = 0; phase < np; ++phase) { - props_[phase] ->rvSat(n, &p[0], &rv[0], &drvdp[0]); + props_[phase] ->rvSat(n, &pvtTableIdx[0], &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]); @@ -262,6 +252,7 @@ BOOST_AUTO_TEST_CASE(test_liveoil) // approximation of the derivatives. const int n = 6; const int np = phase_usage_.num_phases; + std::vector pvtRegionIdx(n, 0); // the tolerance for acceptable difference in values const double reltol = 1e-9; @@ -306,13 +297,13 @@ BOOST_AUTO_TEST_CASE(test_liveoil) } - testmu(reltol, n, np, p, r,z, props_, condition); + testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition); - testb(reltol,n,np,p,r,z,props_,condition); + testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition); - testrsSat(reltol,n,np,p,props_); + testrsSat(reltol,n,np,pvtRegionIdx,p,props_); - testrvSat(reltol,n,np,p,props_); + testrvSat(reltol,n,np,pvtRegionIdx,p,props_); } @@ -337,6 +328,7 @@ BOOST_AUTO_TEST_CASE(test_wetgas) // approximation of the derivatives. const int n = 6; const int np = phase_usage_.num_phases; + std::vector pvtRegionIdx(n, 0); // the tolerance for acceptable difference in values const double reltol = 1e-9; @@ -381,12 +373,12 @@ BOOST_AUTO_TEST_CASE(test_wetgas) } - testmu(reltol, n, np, p, r,z, props_, condition); + testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition); - testb(reltol,n,np,p,r,z,props_,condition); + testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition); - testrsSat(reltol,n,np,p,props_); + testrsSat(reltol,n,np,pvtRegionIdx,p,props_); - testrvSat(reltol,n,np,p,props_); + testrvSat(reltol,n,np,pvtRegionIdx,p,props_); }