diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f86e19cc..3e5db299 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -84,11 +84,11 @@ list (APPEND MAIN_SOURCE_FILES opm/core/props/pvt/BlackoilPvtProperties.cpp opm/core/props/pvt/PvtPropertiesBasic.cpp opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp - opm/core/props/pvt/SinglePvtDead.cpp - opm/core/props/pvt/SinglePvtDeadSpline.cpp - opm/core/props/pvt/SinglePvtInterface.cpp - opm/core/props/pvt/SinglePvtLiveGas.cpp - opm/core/props/pvt/SinglePvtLiveOil.cpp + opm/core/props/pvt/PvtDead.cpp + opm/core/props/pvt/PvtDeadSpline.cpp + opm/core/props/pvt/PvtInterface.cpp + opm/core/props/pvt/PvtLiveGas.cpp + opm/core/props/pvt/PvtLiveOil.cpp opm/core/props/rock/RockBasic.cpp opm/core/props/rock/RockCompressibility.cpp opm/core/props/rock/RockFromDeck.cpp @@ -311,12 +311,12 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/props/pvt/BlackoilPvtProperties.hpp opm/core/props/pvt/PvtPropertiesBasic.hpp opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp - opm/core/props/pvt/SinglePvtConstCompr.hpp - opm/core/props/pvt/SinglePvtDead.hpp - opm/core/props/pvt/SinglePvtDeadSpline.hpp - opm/core/props/pvt/SinglePvtInterface.hpp - opm/core/props/pvt/SinglePvtLiveGas.hpp - opm/core/props/pvt/SinglePvtLiveOil.hpp + opm/core/props/pvt/PvtConstCompr.hpp + opm/core/props/pvt/PvtDead.hpp + opm/core/props/pvt/PvtDeadSpline.hpp + opm/core/props/pvt/PvtInterface.hpp + opm/core/props/pvt/PvtLiveGas.hpp + opm/core/props/pvt/PvtLiveOil.hpp opm/core/props/rock/RockBasic.hpp opm/core/props/rock/RockCompressibility.hpp opm/core/props/rock/RockFromDeck.hpp diff --git a/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp b/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp index 9a9b8bad..713ede5e 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp @@ -23,8 +23,14 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock) { - init(deck, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension, - param, init_rock); + init(deck, + number_of_cells, + global_cell, + cart_dims, + begin_cell_centroids, + dimension, + param, + init_rock); } template @@ -39,7 +45,7 @@ namespace Opm if (init_rock){ rock_.init(deck, number_of_cells, global_cell, cart_dims); } - pvt_.init(deck, /*numSamples=*/0); + pvt_.init(deck, /*numSamples=*/0, number_of_cells, global_cell); SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); satprops_.reset(ptr); @@ -68,7 +74,7 @@ namespace Opm } const int pvt_samples = param.getDefault("pvt_tab_size", 200); - pvt_.init(deck, pvt_samples); + pvt_.init(deck, pvt_samples, number_of_cells, global_cell); // Unfortunate lack of pointer smartness here... const int sat_samples = param.getDefault("sat_tab_size", 200); diff --git a/opm/core/props/pvt/BlackoilPvtProperties.cpp b/opm/core/props/pvt/BlackoilPvtProperties.cpp index 04e56ad1..cb8121a6 100644 --- a/opm/core/props/pvt/BlackoilPvtProperties.cpp +++ b/opm/core/props/pvt/BlackoilPvtProperties.cpp @@ -21,11 +21,11 @@ #include "config.h" #include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include #include #include #include @@ -42,93 +42,107 @@ namespace Opm { } - void BlackoilPvtProperties::init(Opm::DeckConstPtr deck, int samples) + void BlackoilPvtProperties::init(Opm::DeckConstPtr deck, + int numSamples, + int numCompressedCells, + const int *compressedToCartesianCellIdx) { - // If we need multiple regions, this class and the SinglePvt* classes must change. - region_number_ = 0; - phase_usage_ = phaseUsageFromDeck(deck); // Surface densities. Accounting for different orders in eclipse and our code. - if (deck->hasKeyword("DENSITY")) { - Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY"); + Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY"); + int numRegions = densityKeyword->size(); + + densities_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { if (phase_usage_.phase_used[Liquid]) { - densities_[phase_usage_.phase_pos[Liquid]] - = densityKeyword->getRecord(region_number_)->getItem("OIL")->getSIDouble(0); + densities_[regionIdx][phase_usage_.phase_pos[Liquid]] + = densityKeyword->getRecord(regionIdx)->getItem("OIL")->getSIDouble(0); } if (phase_usage_.phase_used[Aqua]) { - densities_[phase_usage_.phase_pos[Aqua]] - = densityKeyword->getRecord(region_number_)->getItem("WATER")->getSIDouble(0); + densities_[regionIdx][phase_usage_.phase_pos[Aqua]] + = densityKeyword->getRecord(regionIdx)->getItem("WATER")->getSIDouble(0); } if (phase_usage_.phase_used[Vapour]) { - densities_[phase_usage_.phase_pos[Vapour]] - = densityKeyword->getRecord(region_number_)->getItem("GAS")->getSIDouble(0); + densities_[regionIdx][phase_usage_.phase_pos[Vapour]] + = densityKeyword->getRecord(regionIdx)->getItem("GAS")->getSIDouble(0); } - } else { - OPM_THROW(std::runtime_error, "Input is missing DENSITY\n"); } - // Set the properties. + // 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); + // Water PVT if (phase_usage_.phase_used[Aqua]) { - if (deck->hasKeyword("PVTW")) { - Opm::PvtwTable pvtwTable(deck->getKeyword("PVTW")); + // 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_); - props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable)); - } else { - // Eclipse 100 default. - props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise)); - } + props_[phase_usage_.phase_pos[Aqua]] = pvtw; } // Oil PVT if (phase_usage_.phase_used[Liquid]) { + // for oil, we support the "PVDO", "PVTO" and "PVCDO" + // keywords... if (deck->hasKeyword("PVDO")) { - Opm::PvdoTable pvdoTable(deck->getKeyword("PVDO"), region_number_); - if (samples > 0) { - props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples)); + Opm::DeckKeywordConstPtr pvdoKeyword = deck->getKeyword("PVDO"); + if (numSamples > 0) { + auto splinePvt = std::shared_ptr(new PvtDeadSpline); + splinePvt->initFromOil(pvdoKeyword, pvtTableIdx_, numSamples); + props_[phase_usage_.phase_pos[Liquid]] = splinePvt; } else { - props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable)); + auto deadPvt = std::shared_ptr(new PvtDead); + deadPvt->initFromOil(pvdoKeyword, pvtTableIdx_); + props_[phase_usage_.phase_pos[Liquid]] = deadPvt; } } else if (deck->hasKeyword("PVTO")) { - Opm::PvtoTable pvtoTable(deck->getKeyword("PVTO"), /*tableIdx=*/0); - - props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable)); + props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"), pvtTableIdx_)); } else if (deck->hasKeyword("PVCDO")) { - Opm::PvcdoTable pvcdoTable(deck->getKeyword("PVCDO")); + std::shared_ptr pvcdo(new PvtConstCompr); + pvcdo->initFromOil(deck->getKeyword("PVCDO"), pvtTableIdx_); - props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable)); + props_[phase_usage_.phase_pos[Liquid]] = pvcdo; } else { - OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n"); + OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n"); } } // Gas PVT if (phase_usage_.phase_used[Vapour]) { + // gas can be specified using the "PVDG" or "PVTG" keywords... if (deck->hasKeyword("PVDG")) { - Opm::PvdgTable pvdgTable(deck->getKeyword("PVDG"), region_number_); - if (samples > 0) { - props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples)); + Opm::DeckKeywordConstPtr pvdgKeyword = deck->getKeyword("PVDG"); + + if (numSamples > 0) { + std::shared_ptr splinePvt(new PvtDeadSpline); + splinePvt->initFromGas(pvdgKeyword, pvtTableIdx_, numSamples); + + props_[phase_usage_.phase_pos[Vapour]] = splinePvt; } else { - props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable)); + std::shared_ptr deadPvt(new PvtDead); + deadPvt->initFromGas(pvdgKeyword, pvtTableIdx_); + + props_[phase_usage_.phase_pos[Vapour]] = deadPvt; } } else if (deck->hasKeyword("PVTG")) { - Opm::PvtgTable pvtgTable(deck->getKeyword("PVTG"), /*tableIdx=*/0); - - props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable)); + props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG"), pvtTableIdx_)); } else { OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n"); } } - - // Must inform pvt property objects of phase structure. - for (int i = 0; i < phase_usage_.num_phases; ++i) { - props_[i]->setPhaseConfiguration(phase_usage_.num_phases, phase_usage_.phase_pos); - } } - const double* BlackoilPvtProperties::surfaceDensities() const + const double* BlackoilPvtProperties::surfaceDensities(int regionIdx) const { - return densities_; + return &densities_[regionIdx][0]; } @@ -234,5 +248,4 @@ namespace Opm } } } - } // namespace Opm diff --git a/opm/core/props/pvt/BlackoilPvtProperties.hpp b/opm/core/props/pvt/BlackoilPvtProperties.hpp index b8bff837..82921c50 100644 --- a/opm/core/props/pvt/BlackoilPvtProperties.hpp +++ b/opm/core/props/pvt/BlackoilPvtProperties.hpp @@ -20,10 +20,11 @@ #ifndef OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED #define OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED -#include +#include #include #include +#include #include #include @@ -37,9 +38,9 @@ namespace Opm /// 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 SinglePvtInterface is that this collects all phases' properties, + /// by PvtInterface is that this collects all phases' properties, /// and therefore the output arrays are of size n*num_phases as opposed - /// to size n in SinglePvtInterface. + /// to size n in PvtInterface. class BlackoilPvtProperties : public BlackoilPhases { public: @@ -47,8 +48,17 @@ namespace Opm BlackoilPvtProperties(); /// 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); + void init(Opm::DeckConstPtr deck, + int samples, + int numCompressedCells, + const int *compressedToCartesianCellIdx); /// \return Object describing the active phases. PhaseUsage phaseUsage() const; @@ -68,7 +78,7 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of size numPhases(). - const double* surfaceDensities() const; + const double* surfaceDensities(int regionIdx = 0) const; /// Viscosity as a function of p and z. void mu(const int n, @@ -109,11 +119,19 @@ namespace Opm PhaseUsage phase_usage_; - int region_number_; + // The PVT properties. One object per active fluid phase. + std::vector > props_; - 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_; - double densities_[MaxNumPhases]; + 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 new file mode 100644 index 00000000..e7606a06 --- /dev/null +++ b/opm/core/props/pvt/PvtConstCompr.hpp @@ -0,0 +1,291 @@ +/* + Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_PVTCONSTCOMPR_HEADER_INCLUDED +#define OPM_PVTCONSTCOMPR_HEADER_INCLUDED + +#include +#include + +#include +#include + +#include +#include + + +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. + class PvtConstCompr : public PvtInterface + { + public: + PvtConstCompr() + {} + + void initFromWater(Opm::DeckKeywordConstPtr pvtwKeyword, const std::vector &pvtTableIdx) + { + pvtTableIdx_ = pvtTableIdx; + + int numRegions = pvtwKeyword->size(); + + ref_press_.resize(numRegions); + ref_B_.resize(numRegions); + comp_.resize(numRegions); + viscosity_.resize(numRegions); + visc_comp_.resize(numRegions); + + for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) { + Opm::DeckRecordConstPtr pvtwRecord = pvtwKeyword->getRecord(regionIdx); + + ref_press_[regionIdx] = pvtwRecord->getItem("P_REF")->getSIDouble(0); + ref_B_[regionIdx] = pvtwRecord->getItem("WATER_VOL_FACTOR")->getSIDouble(0); + comp_[regionIdx] = pvtwRecord->getItem("WATER_COMPRESSIBILITY")->getSIDouble(0); + viscosity_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSITY")->getSIDouble(0); + visc_comp_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSIBILITY")->getSIDouble(0); + } + } + + void initFromOil(Opm::DeckKeywordConstPtr pvcdoKeyword, const std::vector &pvtTableIdx) + { + pvtTableIdx_ = pvtTableIdx; + + int numRegions = pvcdoKeyword->size(); + + ref_press_.resize(numRegions); + ref_B_.resize(numRegions); + comp_.resize(numRegions); + viscosity_.resize(numRegions); + visc_comp_.resize(numRegions); + + for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) { + Opm::DeckRecordConstPtr pvcdoRecord = pvcdoKeyword->getRecord(regionIdx); + + ref_press_[regionIdx] = pvcdoRecord->getItem("P_REF")->getSIDouble(0); + ref_B_[regionIdx] = pvcdoRecord->getItem("BO_REF")->getSIDouble(0); + comp_[regionIdx] = pvcdoRecord->getItem("C_OIL_REF")->getSIDouble(0); + viscosity_[regionIdx] = pvcdoRecord->getItem("MUO_REF")->getSIDouble(0); + visc_comp_[regionIdx] = pvcdoRecord->getItem("OIL_VISCOSIBILITY")->getSIDouble(0); + } + } + + PvtConstCompr(double visc) + : ref_press_(1, 0.0), + ref_B_(1, 1.0), + comp_(1, 0.0), + viscosity_(1, visc), + visc_comp_(1, 0.0) + { + } + + virtual ~PvtConstCompr() + { + } + + virtual void mu(const int n, + const double* p, + const double* /*z*/, + double* output_mu) const + { +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + // Computing a polynomial approximation to the exponential. + int tableIdx = getTableIndex_(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 double* p, + const double* /*r*/, + double* output_mu, + double* output_dmudp, + double* output_dmudr) const + { + // #pragma omp parallel for + for (int i = 0; i < n; ++i) { + // Computing a polynomial approximation to the exponential. + int tableIdx = getTableIndex_(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; + output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx]; + } + std::fill(output_dmudr, output_dmudr + n, 0.0); + } + + virtual void mu(const int n, + const double* p, + const double* /*r*/, + const PhasePresence* /*cond*/, + double* output_mu, + double* output_dmudp, + double* output_dmudr) const + { + // #pragma omp parallel for + for (int i = 0; i < n; ++i) { + // Computing a polynomial approximation to the exponential. + int tableIdx = getTableIndex_(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; + output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx]; + } + std::fill(output_dmudr, output_dmudr + n, 0.0); + } + + virtual void B(const int n, + const double* p, + const double* /*z*/, + double* output_B) const + { +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + // Computing a polynomial approximation to the exponential. + int tableIdx = getTableIndex_(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 double* p, + const double* /*z*/, + double* output_B, + double* output_dBdp) const + { +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + int tableIdx = getTableIndex_(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; + output_dBdp[i] = (-ref_B_[tableIdx]/(d*d))*(1 + x) * comp_[tableIdx]; + } + } + + virtual void b(const int n, + const double* p, + const double* /*r*/, + double* output_b, + double* output_dbdp, + double* output_dbdr) const + { + // #pragma omp parallel for + for (int i = 0; i < n; ++i) { + // Computing a polynomial approximation to the exponential. + int tableIdx = getTableIndex_(i); + double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); + double d = (1.0 + x + 0.5*x*x); + + // b = 1/B = d/ref_B_B; + output_b[i] = d/ref_B_[tableIdx]; + output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx]; + } + + std::fill(output_dbdr, output_dbdr + n, 0.0); + } + + virtual void b(const int n, + const double* p, + const double* /*r*/, + const PhasePresence* /*cond*/, + double* output_b, + double* output_dbdp, + double* output_dbdr) const + { + // #pragma omp parallel for + for (int i = 0; i < n; ++i) { + // Computing a polynomial approximation to the exponential. + int tableIdx = getTableIndex_(i); + double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); + double d = (1.0 + x + 0.5*x*x); + + // b = 1/B = d/ref_B_[tableIdx]B; + output_b[i] = d/ref_B_[tableIdx]; + output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx]; + } + std::fill(output_dbdr, output_dbdr + n, 0.0); + } + + virtual void rsSat(const int n, + const double* /*p*/, + double* output_rsSat, + double* output_drsSatdp) const + { + std::fill(output_rsSat, output_rsSat + n, 0.0); + std::fill(output_drsSatdp, output_drsSatdp + n, 0.0); + } + + virtual void rvSat(const int n, + const double* /*p*/, + double* output_rvSat, + double* output_drvSatdp) const + { + std::fill(output_rvSat, output_rvSat + n, 0.0); + std::fill(output_drvSatdp, output_drvSatdp + n, 0.0); + } + + virtual void R(const int n, + const double* /*p*/, + const double* /*z*/, + double* output_R) const + { + std::fill(output_R, output_R + n, 0.0); + } + + virtual void dRdp(const int n, + const double* /*p*/, + const double* /*z*/, + double* output_R, + double* output_dRdp) const + { + std::fill(output_R, output_R + n, 0.0); + std::fill(output_dRdp, output_dRdp + n, 0.0); + } + + private: + int getTableIndex_(int cellIdx) const + { + if (pvtTableIdx_.empty()) + return 0; + return pvtTableIdx_[cellIdx]; + } + + std::vector pvtTableIdx_; + std::vector ref_press_; + std::vector ref_B_; + std::vector comp_; + std::vector viscosity_; + std::vector visc_comp_; + }; + +} + + +#endif // OPM_PVTCONSTCOMPR_HEADER_INCLUDED + diff --git a/opm/core/props/pvt/SinglePvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp similarity index 58% rename from opm/core/props/pvt/SinglePvtDead.cpp rename to opm/core/props/pvt/PvtDead.cpp index 68449de2..ecd4e885 100644 --- a/opm/core/props/pvt/SinglePvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -19,7 +19,7 @@ #include "config.h" -#include +#include #include // Extra includes for debug dumping of tables. @@ -34,72 +34,84 @@ namespace Opm // Member functions //------------------------------------------------------------------------- /// Constructor - SinglePvtDead::SinglePvtDead(const Opm::PvdoTable& pvdoTable) + void PvtDead::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword, + const std::vector &pvtTableIdx) { - // Copy data - const std::vector& press = pvdoTable.getPressureColumn(); - const std::vector& b = pvdoTable.getFormationFactorColumn(); - const std::vector& visc = pvdoTable.getViscosityColumn(); + pvtTableIdx_ = pvtTableIdx; - const int sz = b.size(); - std::vector bInv(sz); - for (int i = 0; i < sz; ++i) { - bInv[i] = 1.0 / b[i]; + int numRegions = Opm::PvdoTable::numTables(pvdoKeyword); + + // resize the attributes of the object + b_.resize(numRegions); + viscosity_.resize(numRegions); + + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx); + // Copy data + const std::vector& press = pvdoTable.getPressureColumn(); + const std::vector& b = pvdoTable.getFormationFactorColumn(); + const std::vector& visc = pvdoTable.getViscosityColumn(); + + const int sz = b.size(); + std::vector bInv(sz); + for (int i = 0; i < sz; ++i) { + bInv[i] = 1.0 / b[i]; + } + b_[regionIdx] = NonuniformTableLinear(press, bInv); + viscosity_[regionIdx] = NonuniformTableLinear(press, visc); } - b_ = NonuniformTableLinear(press, bInv); - viscosity_ = NonuniformTableLinear(press, visc); - - // Dumping the created tables. -// static int count = 0; -// std::ofstream os((std::string("dump-") + boost::lexical_cast(count++)).c_str()); -// os.precision(15); -// os << "1/B\n\n" << one_over_B_ -// << "\n\nvisc\n\n" << viscosity_ << std::endl; } - /// Constructor - SinglePvtDead::SinglePvtDead(const Opm::PvdgTable& pvdgTable) + + void PvtDead::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword, + const std::vector &pvtTableIdx) { - // Copy data - const std::vector& press = pvdgTable.getPressureColumn(); - const std::vector& b = pvdgTable.getFormationFactorColumn(); - const std::vector& visc = pvdgTable.getViscosityColumn(); + pvtTableIdx_ = pvtTableIdx; - const int sz = b.size(); - std::vector bInv(sz); - for (int i = 0; i < sz; ++i) { - bInv[i] = 1.0 / b[i]; + int numRegions = Opm::PvdgTable::numTables(pvdgKeyword); + + // resize the attributes of the object + b_.resize(numRegions); + viscosity_.resize(numRegions); + + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx); + + // Copy data + const std::vector& press = pvdgTable.getPressureColumn(); + const std::vector& b = pvdgTable.getFormationFactorColumn(); + const std::vector& visc = pvdgTable.getViscosityColumn(); + + const int sz = b.size(); + std::vector bInv(sz); + for (int i = 0; i < sz; ++i) { + bInv[i] = 1.0 / b[i]; + } + b_[regionIdx] = NonuniformTableLinear(press, bInv); + viscosity_[regionIdx] = NonuniformTableLinear(press, visc); } - b_ = NonuniformTableLinear(press, bInv); - viscosity_ = NonuniformTableLinear(press, visc); - - // Dumping the created tables. -// static int count = 0; -// std::ofstream os((std::string("dump-") + boost::lexical_cast(count++)).c_str()); -// os.precision(15); -// os << "1/B\n\n" << one_over_B_ -// << "\n\nvisc\n\n" << viscosity_ << std::endl; } // Destructor - SinglePvtDead::~SinglePvtDead() + PvtDead::~PvtDead() { } - void SinglePvtDead::mu(const int n, + void PvtDead::mu(const int n, const double* p, const double* /*z*/, double* output_mu) const { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - output_mu[i] = viscosity_(p[i]); + int regionIdx = pvtTableIdx_[i]; + output_mu[i] = viscosity_[regionIdx](p[i]); } } - void SinglePvtDead::mu(const int n, + void PvtDead::mu(const int n, const double* p, const double* /*r*/, double* output_mu, @@ -108,14 +120,15 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - output_mu[i] = viscosity_(p[i]); - output_dmudp[i] = viscosity_.derivative(p[i]); + int regionIdx = pvtTableIdx_[i]; + output_mu[i] = viscosity_[regionIdx](p[i]); + output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); } std::fill(output_dmudr, output_dmudr + n, 0.0); } - void SinglePvtDead::mu(const int n, + void PvtDead::mu(const int n, const double* p, const double* /*r*/, const PhasePresence* /*cond*/, @@ -125,14 +138,15 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - output_mu[i] = viscosity_(p[i]); - output_dmudp[i] = viscosity_.derivative(p[i]); + int regionIdx = pvtTableIdx_[i]; + output_mu[i] = viscosity_[regionIdx](p[i]); + output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); } std::fill(output_dmudr, output_dmudr + n, 0.0); } - void SinglePvtDead::B(const int n, + void PvtDead::B(const int n, const double* p, const double* /*z*/, double* output_B) const @@ -140,11 +154,12 @@ namespace Opm // #pragma omp parallel for // B = 1/b for (int i = 0; i < n; ++i) { - output_B[i] = 1.0/b_(p[i]); + int regionIdx = pvtTableIdx_[i]; + output_B[i] = 1.0/b_[regionIdx](p[i]); } } - void SinglePvtDead::dBdp(const int n, + void PvtDead::dBdp(const int n, const double* p, const double* /*z*/, double* output_B, @@ -153,12 +168,13 @@ namespace Opm B(n, p, 0, output_B); // #pragma omp parallel for for (int i = 0; i < n; ++i) { + int regionIdx = pvtTableIdx_[i]; double Bg = output_B[i]; - output_dBdp[i] = -Bg*Bg*b_.derivative(p[i]); + output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]); } } - void SinglePvtDead::b(const int n, + void PvtDead::b(const int n, const double* p, const double* /*r*/, double* output_b, @@ -168,15 +184,17 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - output_b[i] = b_(p[i]); - output_dbdp[i] = b_.derivative(p[i]); + int regionIdx = pvtTableIdx_[i]; + + output_b[i] = b_[regionIdx](p[i]); + output_dbdp[i] = b_[regionIdx].derivative(p[i]); } std::fill(output_dbdr, output_dbdr + n, 0.0); } - void SinglePvtDead::b(const int n, + void PvtDead::b(const int n, const double* p, const double* /*r*/, const PhasePresence* /*cond*/, @@ -187,15 +205,17 @@ namespace Opm { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - output_b[i] = b_(p[i]); - output_dbdp[i] = b_.derivative(p[i]); + int regionIdx = pvtTableIdx_[i]; + + output_b[i] = b_[regionIdx](p[i]); + output_dbdp[i] = b_[regionIdx].derivative(p[i]); } std::fill(output_dbdr, output_dbdr + n, 0.0); } - void SinglePvtDead::rsSat(const int n, + void PvtDead::rsSat(const int n, const double* /*p*/, double* output_rsSat, double* output_drsSatdp) const @@ -204,7 +224,7 @@ namespace Opm std::fill(output_drsSatdp, output_drsSatdp + n, 0.0); } - void SinglePvtDead::rvSat(const int n, + void PvtDead::rvSat(const int n, const double* /*p*/, double* output_rvSat, double* output_drvSatdp) const @@ -213,7 +233,7 @@ namespace Opm std::fill(output_drvSatdp, output_drvSatdp + n, 0.0); } - void SinglePvtDead::R(const int n, + void PvtDead::R(const int n, const double* /*p*/, const double* /*z*/, double* output_R) const @@ -221,7 +241,7 @@ namespace Opm std::fill(output_R, output_R + n, 0.0); } - void SinglePvtDead::dRdp(const int n, + void PvtDead::dRdp(const int n, const double* /*p*/, const double* /*z*/, double* output_R, diff --git a/opm/core/props/pvt/SinglePvtDead.hpp b/opm/core/props/pvt/PvtDead.hpp similarity index 88% rename from opm/core/props/pvt/SinglePvtDead.hpp rename to opm/core/props/pvt/PvtDead.hpp index e7d28872..c7d26e6b 100644 --- a/opm/core/props/pvt/SinglePvtDead.hpp +++ b/opm/core/props/pvt/PvtDead.hpp @@ -17,11 +17,11 @@ along with OPM. If not, see . */ -#ifndef OPM_SINGLEPVTDEAD_HEADER_INCLUDED -#define OPM_SINGLEPVTDEAD_HEADER_INCLUDED +#ifndef OPM_PVTDEAD_HEADER_INCLUDED +#define OPM_PVTDEAD_HEADER_INCLUDED -#include +#include #include #include @@ -39,12 +39,16 @@ namespace Opm /// 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. - class SinglePvtDead : public SinglePvtInterface + class PvtDead : public PvtInterface { public: - SinglePvtDead(const Opm::PvdoTable& pvdoTable); - SinglePvtDead(const Opm::PvdgTable& pvdgTable); - virtual ~SinglePvtDead(); + PvtDead() {}; + + void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword, + const std::vector &pvtTableIdx); + void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword, + const std::vector &pvtTableIdx); + virtual ~PvtDead(); /// Viscosity as a function of p and z. virtual void mu(const int n, @@ -130,11 +134,11 @@ namespace Opm double* output_dRdp) const; private: // PVT properties of dry gas or dead oil - NonuniformTableLinear b_; - NonuniformTableLinear viscosity_; + std::vector pvtTableIdx_; + std::vector > b_; + std::vector > viscosity_; }; - } -#endif // OPM_SINGLEPVTDEAD_HEADER_INCLUDED +#endif // OPM_PVTDEAD_HEADER_INCLUDED diff --git a/opm/core/props/pvt/PvtDeadSpline.cpp b/opm/core/props/pvt/PvtDeadSpline.cpp new file mode 100644 index 00000000..1b40fae4 --- /dev/null +++ b/opm/core/props/pvt/PvtDeadSpline.cpp @@ -0,0 +1,255 @@ +/* + Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include "config.h" +#include +#include + +#include + +#include + +// Extra includes for debug dumping of tables. +// #include +// #include +// #include + +namespace Opm +{ + //------------------------------------------------------------------------ + // Member functions + //------------------------------------------------------------------------- + + PvtDeadSpline::PvtDeadSpline() + {} + + 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 + b_.resize(numRegions); + viscosity_.resize(numRegions); + + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx); + + int numRows = pvdoTable.numRows(); + + // Copy data + const std::vector &press = pvdoTable.getPressureColumn(); + const std::vector &B = pvdoTable.getFormationFactorColumn(); + const std::vector &visc = pvdoTable.getViscosityColumn(); + + std::vector B_inv(numRows); + for (int i = 0; i < numRows; ++i) { + B_inv[i] = 1.0 / B[i]; + } + + buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]); + buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]); + } + } + + 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 + b_.resize(numRegions); + viscosity_.resize(numRegions); + + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx); + + int numRows = pvdgTable.numRows(); + + // Copy data + const std::vector &press = pvdgTable.getPressureColumn(); + const std::vector &B = pvdgTable.getFormationFactorColumn(); + const std::vector &visc = pvdgTable.getViscosityColumn(); + + std::vector B_inv(numRows); + for (int i = 0; i < numRows; ++i) { + B_inv[i] = 1.0 / B[i]; + } + + buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]); + buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]); + } + } + + // Destructor + PvtDeadSpline::~PvtDeadSpline() + { + } + + + + void PvtDeadSpline::mu(const int n, + 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]; + output_mu[i] = viscosity_[regionIdx](p[i]); + } + } + + void PvtDeadSpline::mu(const int n, + const double* p, + const double* /*r*/, + double* output_mu, + double* output_dmudp, + double* output_dmudr) const + { + // #pragma omp parallel for + for (int i = 0; i < n; ++i) { + int regionIdx = pvtTableIdx_[i]; + output_mu[i] = viscosity_[regionIdx](p[i]); + output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); + } + std::fill(output_dmudr, output_dmudr + n, 0.0); + } + + void PvtDeadSpline::mu(const int n, + const double* p, + const double* /*r*/, + const PhasePresence* /*cond*/, + double* output_mu, + double* output_dmudp, + double* output_dmudr) const + { + // #pragma omp parallel for + + for (int i = 0; i < n; ++i) { + int regionIdx = pvtTableIdx_[i]; + output_mu[i] = viscosity_[regionIdx](p[i]); + output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); + } + std::fill(output_dmudr, output_dmudr + n, 0.0); + } + + void PvtDeadSpline::B(const int n, + 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]; + output_B[i] = 1.0/b_[regionIdx](p[i]); + } + } + + void PvtDeadSpline::dBdp(const int n, + const double* p, + const double* /*z*/, + double* output_B, + double* output_dBdp) const + { + B(n, p, 0, output_B); +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + int regionIdx = pvtTableIdx_[i]; + double Bg = output_B[i]; + output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]); + } + } + + void PvtDeadSpline::b(const int n, + const double* p, + const double* /*r*/, + double* output_b, + double* output_dbdp, + double* output_dbdr) const + { + // #pragma omp parallel for + for (int i = 0; i < n; ++i) { + int regionIdx = pvtTableIdx_[i]; + output_b[i] = b_[regionIdx](p[i]); + output_dbdp[i] = b_[regionIdx].derivative(p[i]); + } + std::fill(output_dbdr, output_dbdr + n, 0.0); + } + + void PvtDeadSpline::b(const int n, + const double* p, + const double* /*r*/, + const PhasePresence* /*cond*/, + double* output_b, + double* output_dbdp, + double* output_dbdr) const + { + // #pragma omp parallel for + for (int i = 0; i < n; ++i) { + int regionIdx = pvtTableIdx_[i]; + output_b[i] = b_[regionIdx](p[i]); + output_dbdp[i] = b_[regionIdx].derivative(p[i]); + } + std::fill(output_dbdr, output_dbdr + n, 0.0); + } + + void PvtDeadSpline::rsSat(const int n, + const double* /*p*/, + double* output_rsSat, + double* output_drsSatdp) const + { + std::fill(output_rsSat, output_rsSat + n, 0.0); + std::fill(output_drsSatdp, output_drsSatdp + n, 0.0); + } + + void PvtDeadSpline::rvSat(const int n, + const double* /*p*/, + double* output_rvSat, + double* output_drvSatdp) const + { + std::fill(output_rvSat, output_rvSat + n, 0.0); + std::fill(output_drvSatdp, output_drvSatdp + n, 0.0); + } + + void PvtDeadSpline::R(const int n, + const double* /*p*/, + const double* /*z*/, + double* output_R) const + { + std::fill(output_R, output_R + n, 0.0); + } + + void PvtDeadSpline::dRdp(const int n, + const double* /*p*/, + const double* /*z*/, + double* output_R, + double* output_dRdp) const + { + std::fill(output_R, output_R + n, 0.0); + std::fill(output_dRdp, output_dRdp + n, 0.0); + } + +} diff --git a/opm/core/props/pvt/SinglePvtDeadSpline.hpp b/opm/core/props/pvt/PvtDeadSpline.hpp similarity index 86% rename from opm/core/props/pvt/SinglePvtDeadSpline.hpp rename to opm/core/props/pvt/PvtDeadSpline.hpp index 88320ed2..b47a4986 100644 --- a/opm/core/props/pvt/SinglePvtDeadSpline.hpp +++ b/opm/core/props/pvt/PvtDeadSpline.hpp @@ -17,10 +17,10 @@ along with OPM. If not, see . */ -#ifndef OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED -#define OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED +#ifndef OPM_PVTDEADSPLINE_HEADER_INCLUDED +#define OPM_PVTDEADSPLINE_HEADER_INCLUDED -#include +#include #include #include @@ -38,12 +38,19 @@ namespace Opm /// 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. - class SinglePvtDeadSpline : public SinglePvtInterface + class PvtDeadSpline : public PvtInterface { public: - SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples); - SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples); - virtual ~SinglePvtDeadSpline(); + 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, @@ -128,11 +135,12 @@ namespace Opm double* output_dRdp) const; private: // PVT properties of dry gas or dead oil - UniformTableLinear b_; - UniformTableLinear viscosity_; + std::vector pvtTableIdx_; + std::vector > b_; + std::vector > viscosity_; }; } -#endif // OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED +#endif // OPM_PVTDEADSPLINE_HEADER_INCLUDED diff --git a/opm/core/props/pvt/PvtInterface.cpp b/opm/core/props/pvt/PvtInterface.cpp new file mode 100644 index 00000000..94852f4f --- /dev/null +++ b/opm/core/props/pvt/PvtInterface.cpp @@ -0,0 +1,79 @@ +/* + Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include "config.h" +#include + +#include + +namespace Opm +{ + + PvtInterface::PvtInterface() + : num_phases_(MaxNumPhases) + { + for (int i = 0; i < MaxNumPhases; ++i) { + phase_pos_[i] = i; + } + } + + PvtInterface::~PvtInterface() + { + } + + void PvtInterface::setPhaseConfiguration(const int num_phases, const int* phase_pos) + { + num_phases_ = num_phases; + for (int i = 0; i < MaxNumPhases; ++i) { + phase_pos_[i] = phase_pos[i]; + } + } + + void extractPvtTableIndex(std::vector &pvtTableIdx, + Opm::DeckConstPtr deck, + size_t numCompressed, + const int *compressedToCartesianCellIdx) + { + if (deck->hasKeyword("PVTNUM")) { + // we have a "real" multi-pvt deck. + + // get the PVT number from the deck + Opm::DeckKeywordConstPtr pvtnumKeyword = deck->getKeyword("PVTNUM"); + const std::vector &pvtnumData = pvtnumKeyword->getIntData(); + + // convert this into an array of compressed cells. since + // Eclipse uses Fortran-style indices which start at 1 + // instead of 0, we subtract 1. + pvtTableIdx.resize(numCompressed); + for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) { + size_t cartesianCellIdx = compressedToCartesianCellIdx[cellIdx]; + assert(cartesianCellIdx < pvtnumData.size()); + + pvtTableIdx[cellIdx] = pvtnumData[cartesianCellIdx] - 1; + } + } + else { + // create the pvtIdxArray: all cells use the one and only + // PVT table... + pvtTableIdx.resize(numCompressed); + std::fill(pvtTableIdx.begin(), pvtTableIdx.end(), 0); + } + } + +} // namespace Opm diff --git a/opm/core/props/pvt/SinglePvtInterface.hpp b/opm/core/props/pvt/PvtInterface.hpp similarity index 85% rename from opm/core/props/pvt/SinglePvtInterface.hpp rename to opm/core/props/pvt/PvtInterface.hpp index 83d3febb..8be7437a 100644 --- a/opm/core/props/pvt/SinglePvtInterface.hpp +++ b/opm/core/props/pvt/PvtInterface.hpp @@ -17,22 +17,22 @@ along with OPM. If not, see . */ -#ifndef OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED -#define OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED - +#ifndef OPM_PVTINTERFACE_HEADER_INCLUDED +#define OPM_PVTINTERFACE_HEADER_INCLUDED #include +#include namespace Opm { - class SinglePvtInterface : public BlackoilPhases + class PvtInterface : public BlackoilPhases { public: - SinglePvtInterface(); + PvtInterface(); - virtual ~SinglePvtInterface(); + virtual ~PvtInterface(); /// \param[in] num_phases The number of active phases. /// \param[in] phase_pos Array of BlackpoilPhases::MaxNumPhases @@ -138,7 +138,22 @@ namespace Opm int phase_pos_[MaxNumPhases]; }; + /*! + * \brief Helper function to create an array containing the (C-Style) + * PVT table index for each compressed cell. + * + * 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. + */ + void extractPvtTableIndex(std::vector &pvtTableIdx, + Opm::DeckConstPtr deck, + size_t numCompressed, + const int *compressedToCartesianIdx); + } // namespace Opm -#endif // OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED +#endif // OPM_PVTINTERFACE_HEADER_INCLUDED diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp new file mode 100644 index 00000000..6e1109a6 --- /dev/null +++ b/opm/core/props/pvt/PvtLiveGas.cpp @@ -0,0 +1,516 @@ +//=========================================================================== +// +// File: MiscibilityLiveGas.cpp +// +// Created: Wed Feb 10 09:21:53 2010 +// +// Author: Bjørn Spjelkavik +// +// Revision: $Id$ +// +//=========================================================================== +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include "config.h" +#include +#include +#include + +#include +#include + +namespace Opm +{ + + using Opm::linearInterpolation; + using Opm::linearInterpolationDerivative; + + + //------------------------------------------------------------------------ + // Member functions + //------------------------------------------------------------------------- + PvtLiveGas::PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword, + const std::vector &pvtTableIdx) + { + pvtTableIdx_ = pvtTableIdx; + + int numTables = Opm::PvtgTable::numTables(pvtgKeyword); + saturated_gas_table_.resize(numTables); + undersat_gas_tables_.resize(numTables); + + for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) { + Opm::PvtgTable pvtgTable(pvtgKeyword, pvtTableIdx); + + // GAS, PVTG + saturated_gas_table_[pvtTableIdx].resize(4); + saturated_gas_table_[pvtTableIdx][0] = pvtgTable.getOuterTable()->getPressureColumn(); + saturated_gas_table_[pvtTableIdx][1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); + saturated_gas_table_[pvtTableIdx][2] = pvtgTable.getOuterTable()->getGasViscosityColumn(); + saturated_gas_table_[pvtTableIdx][3] = pvtgTable.getOuterTable()->getOilSolubilityColumn(); + + int sz = pvtgTable.getOuterTable()->numRows(); + undersat_gas_tables_[pvtTableIdx].resize(sz); + for (int i=0; i 1/Bg + for (int i=0; i > &saturatedGasTable = + saturated_gas_table_[pvtTableIdx]; + const std::vector > > &undersatGasTables = + undersat_gas_tables_[pvtTableIdx]; + + int section; + double Rval = linearInterpolation(saturatedGasTable[0], + saturatedGasTable[3], press, + section); + double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]]; + if (deriv) { + if (Rval < maxR ) { // Saturated case + return linearInterpolationDerivative(saturatedGasTable[0], + saturatedGasTable[item], + press); + } else { // Undersaturated case + int is = section; + if (undersatGasTables[is][0].size() < 2) { + double val = (saturatedGasTable[item][is+1] + - saturatedGasTable[item][is]) / + (saturatedGasTable[0][is+1] - + saturatedGasTable[0][is]); + return val; + } + double val1 = + linearInterpolation(undersatGasTables[is][0], + undersatGasTables[is][item], + maxR); + double val2 = + linearInterpolation(undersatGasTables[is+1][0], + undersatGasTables[is+1][item], + maxR); + double val = (val2 - val1)/ + (saturatedGasTable[0][is+1] - saturatedGasTable[0][is]); + return val; + } + } else { + if (Rval < maxR ) { // Saturated case + return linearInterpolation(saturatedGasTable[0], + saturatedGasTable[item], + press); + } else { // Undersaturated case + int is = section; + // Extrapolate from first table section + if (is == 0 && press < saturatedGasTable[0][0]) { + return linearInterpolation(undersatGasTables[0][0], + undersatGasTables[0][item], + maxR); + } + + // Extrapolate from last table section + int ltp = saturatedGasTable[0].size() - 1; + if (is+1 == ltp && press > saturatedGasTable[0][ltp]) { + return linearInterpolation(undersatGasTables[ltp][0], + undersatGasTables[ltp][item], + maxR); + } + + // Interpolate between table sections + double w = (press - saturatedGasTable[0][is]) / + (saturatedGasTable[0][is+1] - + saturatedGasTable[0][is]); + if (undersatGasTables[is][0].size() < 2) { + double val = saturatedGasTable[item][is] + + w*(saturatedGasTable[item][is+1] - + saturatedGasTable[item][is]); + return val; + } + double val1 = + linearInterpolation(undersatGasTables[is][0], + undersatGasTables[is][item], + maxR); + double val2 = + linearInterpolation(undersatGasTables[is+1][0], + undersatGasTables[is+1][item], + maxR); + double val = val1 + w*(val2 - val1); + return val; + } + } + } + + double PvtLiveGas::miscible_gas(const double press, + const double r, + const PhasePresence& cond, + const int pvtTableIdx, + const int item, + const int deriv) const + { + const std::vector > &saturatedGasTable = + saturated_gas_table_[pvtTableIdx]; + const std::vector > > &undersatGasTables = + undersat_gas_tables_[pvtTableIdx]; + + const bool isSat = cond.hasFreeOil(); + + // Derivative w.r.t p + if (deriv == 1) { + if (isSat) { // Saturated case + return linearInterpolationDerivative(saturatedGasTable[0], + saturatedGasTable[item], + press); + } else { // Undersaturated case + int is = tableIndex(saturatedGasTable[0], press); + if (undersatGasTables[is][0].size() < 2) { + double val = (saturatedGasTable[item][is+1] + - saturatedGasTable[item][is]) / + (saturatedGasTable[0][is+1] - + saturatedGasTable[0][is]); + return val; + } + double val1 = + linearInterpolation(undersatGasTables[is][0], + undersatGasTables[is][item], + r); + double val2 = + linearInterpolation(undersatGasTables[is+1][0], + undersatGasTables[is+1][item], + r); + double val = (val2 - val1)/ + (saturatedGasTable[0][is+1] - saturatedGasTable[0][is]); + return val; + } + } else if (deriv == 2){ + if (isSat) { + return 0; + } else { + int is = tableIndex(saturatedGasTable[0], press); + double w = (press - saturatedGasTable[0][is]) / + (saturatedGasTable[0][is+1] - saturatedGasTable[0][is]); + assert(undersatGasTables[is][0].size() >= 2); + assert(undersatGasTables[is+1][0].size() >= 2); + double val1 = + linearInterpolationDerivative(undersatGasTables[is][0], + undersatGasTables[is][item], + r); + double val2 = + linearInterpolationDerivative(undersatGasTables[is+1][0], + undersatGasTables[is+1][item], + r); + + double val = val1 + w * (val2 - val1); + return val; + + } + } else { + if (isSat) { // Saturated case + return linearInterpolation(saturatedGasTable[0], + saturatedGasTable[item], + press); + } else { // Undersaturated case + int is = tableIndex(saturatedGasTable[0], press); + // Extrapolate from first table section + if (is == 0 && press < saturatedGasTable[0][0]) { + return linearInterpolation(undersatGasTables[0][0], + undersatGasTables[0][item], + r); + } + + // Extrapolate from last table section + //int ltp = saturatedGasTable[0].size() - 1; + //if (is+1 == ltp && press > saturatedGasTable[0][ltp]) { + // return linearInterpolation(undersatGasTables[ltp][0], + // undersatGasTables[ltp][item], + // r); + //} + + // Interpolate between table sections + double w = (press - saturatedGasTable[0][is]) / + (saturatedGasTable[0][is+1] - + saturatedGasTable[0][is]); + if (undersatGasTables[is][0].size() < 2) { + double val = saturatedGasTable[item][is] + + w*(saturatedGasTable[item][is+1] - + saturatedGasTable[item][is]); + return val; + } + double val1 = + linearInterpolation(undersatGasTables[is][0], + undersatGasTables[is][item], + r); + double val2 = + linearInterpolation(undersatGasTables[is+1][0], + undersatGasTables[is+1][item], + r); + double val = val1 + w*(val2 - val1); + return val; + } + } + } + + + + +} // namespace Opm diff --git a/opm/core/props/pvt/SinglePvtLiveGas.hpp b/opm/core/props/pvt/PvtLiveGas.hpp similarity index 84% rename from opm/core/props/pvt/SinglePvtLiveGas.hpp rename to opm/core/props/pvt/PvtLiveGas.hpp index b61bd3e4..7cdc361e 100644 --- a/opm/core/props/pvt/SinglePvtLiveGas.hpp +++ b/opm/core/props/pvt/PvtLiveGas.hpp @@ -17,10 +17,10 @@ along with OPM. If not, see . */ -#ifndef OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED -#define OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED +#ifndef OPM_PVTLIVEGAS_HEADER_INCLUDED +#define OPM_PVTLIVEGAS_HEADER_INCLUDED -#include +#include #include @@ -35,11 +35,11 @@ namespace Opm /// 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. - class SinglePvtLiveGas : public SinglePvtInterface + class PvtLiveGas : public PvtInterface { public: - SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable); - virtual ~SinglePvtLiveGas(); + PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword, const std::vector &pvtTableIdx); + virtual ~PvtLiveGas(); /// Viscosity as a function of p and z. virtual void mu(const int n, @@ -126,28 +126,30 @@ namespace Opm double* output_dRdp) const; protected: - double evalB(double press, const double* surfvol) const; - void evalBDeriv(double press, const double* surfvol, double& B, double& dBdp) const; - double evalR(double press, const double* surfvol) const; - void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const; + 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; + void evalRDeriv(double press, const double* surfvol, int pvtTableIdx, double& R, double& dRdp) const; // item: 1=>B 2=>mu; double miscible_gas(const double press, const double* surfvol, + const int pvtTableIdx, const int item, const bool deriv = false) const; double miscible_gas(const double press, const double r, const PhasePresence& cond, + const int pvtTableIdx, const int item, const int deriv = 0) const; // PVT properties of wet gas (with vaporised oil) - std::vector > saturated_gas_table_; - std::vector > > undersat_gas_tables_; - + std::vector pvtTableIdx_; + std::vector< std::vector > > saturated_gas_table_; + std::vector< std::vector > > > undersat_gas_tables_; }; } -#endif // OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED +#endif // OPM_PVTLIVEGAS_HEADER_INCLUDED diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp new file mode 100644 index 00000000..2c6691f8 --- /dev/null +++ b/opm/core/props/pvt/PvtLiveOil.cpp @@ -0,0 +1,583 @@ +/* + Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include "config.h" +#include +#include +#include + +#include + +namespace Opm +{ + + using Opm::linearInterpolation; + using Opm::linearInterpolationDerivative; + using Opm::tableIndex; + + //------------------------------------------------------------------------ + // Member functions + //------------------------------------------------------------------------- + PvtLiveOil::PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword, + const std::vector &pvtTableIdx) + { + pvtTableIdx_ = pvtTableIdx; + + int numTables = Opm::PvtoTable::numTables(pvtoKeyword); + saturated_oil_table_.resize(numTables); + undersat_oil_tables_.resize(numTables); + + for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) { + Opm::PvtoTable pvtoTable(pvtoKeyword, pvtTableIdx); + + const auto saturatedPvto = pvtoTable.getOuterTable(); + + // OIL, PVTO + saturated_oil_table_[pvtTableIdx].resize(4); + const int sz = saturatedPvto->numRows(); + for (int k=0; k<4; ++k) { + saturated_oil_table_[pvtTableIdx][k].resize(sz); + } + for (int i=0; igetPressureColumn()[i]; // p + saturated_oil_table_[pvtTableIdx][1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo + saturated_oil_table_[pvtTableIdx][2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o + saturated_oil_table_[pvtTableIdx][3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs + } + + undersat_oil_tables_[pvtTableIdx].resize(sz); + for (int i=0; inumRows(); + undersat_oil_tables_[pvtTableIdx][i][0].resize(tsize); + undersat_oil_tables_[pvtTableIdx][i][1].resize(tsize); + undersat_oil_tables_[pvtTableIdx][i][2].resize(tsize); + for (int j=0; jgetPressureColumn()[j]; // p + undersat_oil_tables_[pvtTableIdx][i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j]; // 1/Bo + undersat_oil_tables_[pvtTableIdx][i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j]; // mu_o + } + } + + // Complete undersaturated tables by extrapolating from existing data + // as is done in Eclipse and Mrst + int iNext = -1; + for (int i=0; i 1) { + continue; + } + // Look ahead for next record containing undersaturated data + if (iNext < i) { + iNext = i+1; + while (iNext > >::size_type sz_t; + for (sz_t j=1; j= 2); + assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); + double val1 = + linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0], + undersat_oil_tables_[pvtTableIdx][is][item], + press); + double val2 = + linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0], + undersat_oil_tables_[pvtTableIdx][is+1][item], + press); + double val = val1 + w*(val2 - val1); + return val; + } + } else { + if (Rval < maxR ) { // Saturated case + return linearInterpolation(saturated_oil_table_[pvtTableIdx][0], + saturated_oil_table_[pvtTableIdx][item], + press); + } else { // Undersaturated case + // Interpolate between table sections + int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], maxR); + double w = (maxR - saturated_oil_table_[pvtTableIdx][3][is]) / + (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]); + assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); + assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); + double val1 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], + undersat_oil_tables_[pvtTableIdx][is][item], + press); + double val2 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], + undersat_oil_tables_[pvtTableIdx][is+1][item], + press); + double val = val1 + w*(val2 - val1); + return val; + } + } + } + + double PvtLiveOil::miscible_oil(const double press, + const double r, + const int pvtTableIdx, + const int item, + const int deriv) const + { + int section; + double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0], + saturated_oil_table_[pvtTableIdx][3], + press, section); + // derivative with respect to frist component (pressure) + if (deriv == 1) { + if (Rval <= r ) { // Saturated case + return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0], + saturated_oil_table_[pvtTableIdx][item], + press); + } else { // Undersaturated case + int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); + double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) / + (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]); + assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); + assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); + double val1 = + linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0], + undersat_oil_tables_[pvtTableIdx][is][item], + press); + double val2 = + linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0], + undersat_oil_tables_[pvtTableIdx][is+1][item], + press); + double val = val1 + w*(val2 - val1); + return val; + } + // derivative with respect to second component (r) + } else if (deriv == 2) { + if (Rval <= r ) { // Saturated case + return 0; + } else { // Undersaturated case + int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); + assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); + assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); + double val1 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], + undersat_oil_tables_[pvtTableIdx][is][item], + press); + double val2 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], + undersat_oil_tables_[pvtTableIdx][is+1][item], + press); + + double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][3][is+1]-saturated_oil_table_[pvtTableIdx][3][is]); + return val; + } + + + } else { + if (Rval <= r ) { // Saturated case + return linearInterpolation(saturated_oil_table_[pvtTableIdx][0], + saturated_oil_table_[pvtTableIdx][item], + press); + } else { // Undersaturated case + // Interpolate between table sections + int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); + double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) / + (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]); + assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); + assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); + double val1 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], + undersat_oil_tables_[pvtTableIdx][is][item], + press); + double val2 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], + undersat_oil_tables_[pvtTableIdx][is+1][item], + press); + double val = val1 + w*(val2 - val1); + return val; + } + } + } + + double PvtLiveOil::miscible_oil(const double press, + const double r, + const PhasePresence& cond, + const int pvtTableIdx, + const int item, + const int deriv) const + { + const bool isSat = cond.hasFreeGas(); + + // derivative with respect to frist component (pressure) + if (deriv == 1) { + if (isSat) { // Saturated case + return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0], + saturated_oil_table_[pvtTableIdx][item], + press); + } else { // Undersaturated case + int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); + double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) / + (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]); + assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); + assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); + double val1 = + linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0], + undersat_oil_tables_[pvtTableIdx][is][item], + press); + double val2 = + linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0], + undersat_oil_tables_[pvtTableIdx][is+1][item], + press); + double val = val1 + w*(val2 - val1); + return val; + } + // derivative with respect to second component (r) + } else if (deriv == 2) { + if (isSat) { // Saturated case + return 0; + } else { // Undersaturated case + int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); + assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); + assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); + double val1 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], + undersat_oil_tables_[pvtTableIdx][is][item], + press); + double val2 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], + undersat_oil_tables_[pvtTableIdx][is+1][item], + press); + + double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][3][is+1]-saturated_oil_table_[pvtTableIdx][3][is]); + return val; + } + + + } else { + if (isSat) { // Saturated case + return linearInterpolation(saturated_oil_table_[pvtTableIdx][0], + saturated_oil_table_[pvtTableIdx][item], + press); + } else { // Undersaturated case + // Interpolate between table sections + int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); + double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) / + (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]); + assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); + assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); + double val1 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], + undersat_oil_tables_[pvtTableIdx][is][item], + press); + double val2 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], + undersat_oil_tables_[pvtTableIdx][is+1][item], + press); + double val = val1 + w*(val2 - val1); + return val; + } + } + } + +} // namespace Opm diff --git a/opm/core/props/pvt/SinglePvtLiveOil.hpp b/opm/core/props/pvt/PvtLiveOil.hpp similarity index 83% rename from opm/core/props/pvt/SinglePvtLiveOil.hpp rename to opm/core/props/pvt/PvtLiveOil.hpp index 8f5328f1..b16708a0 100644 --- a/opm/core/props/pvt/SinglePvtLiveOil.hpp +++ b/opm/core/props/pvt/PvtLiveOil.hpp @@ -17,11 +17,12 @@ along with OPM. If not, see . */ -#ifndef OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED -#define OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED +#ifndef OPM_PVTLIVEOIL_HEADER_INCLUDED +#define OPM_PVTLIVEOIL_HEADER_INCLUDED -#include +#include +#include #include #include @@ -35,11 +36,11 @@ namespace Opm /// 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. - class SinglePvtLiveOil : public SinglePvtInterface + class PvtLiveOil : public PvtInterface { public: - SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable); - virtual ~SinglePvtLiveOil(); + PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword, const std::vector &pvtTableIdx); + virtual ~PvtLiveOil(); /// Viscosity as a function of p and z. virtual void mu(const int n, @@ -124,34 +125,38 @@ namespace Opm double* output_dRdp) const; private: - double evalB(double press, const double* surfvol) const; - void evalBDeriv(double press, const double* surfvol, double& B, double& dBdp) const; - double evalR(double press, const double* surfvol) const; - void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const; + 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; + void evalRDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& R, double& dRdp) const; // item: 1=>1/B 2=>mu; double miscible_oil(const double press, const double* surfvol, + const int pvtTableIdx, const int item, const bool deriv = false) const; double miscible_oil(const double press, const double r, + const int pvtTableIdx, const int item, const int deriv = 0) const; double miscible_oil(const double press, const double r, const PhasePresence& cond, + const int pvtTableIdx, const int item, const int deriv = 0) const; // PVT properties of live oil (with dissolved gas) - std::vector > saturated_oil_table_; - std::vector > > undersat_oil_tables_; + std::vector pvtTableIdx_; + std::vector > > saturated_oil_table_; + std::vector > > > undersat_oil_tables_; }; } -#endif // OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED +#endif // OPM_PVTLIVEOIL_HEADER_INCLUDED diff --git a/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp b/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp index c4c764ac..e23839e9 100644 --- a/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp +++ b/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp @@ -35,7 +35,7 @@ namespace Opm void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck ) { - // If we need multiple regions, this class and the SinglePvt* classes must change. + // If we need multiple regions, this class and the Pvt* classes must change. int region_number = 0; PhaseUsage phase_usage = phaseUsageFromDeck(deck); diff --git a/opm/core/props/pvt/SinglePvtConstCompr.hpp b/opm/core/props/pvt/SinglePvtConstCompr.hpp deleted file mode 100644 index 180611ca..00000000 --- a/opm/core/props/pvt/SinglePvtConstCompr.hpp +++ /dev/null @@ -1,291 +0,0 @@ -/* - Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED -#define OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED - -#include -#include - -#include -#include - -#include -#include - - -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. - class SinglePvtConstCompr : public SinglePvtInterface - { - public: - SinglePvtConstCompr(const Opm::PvtwTable &pvtwTable) - { - if (pvtwTable.numRows() != 1) - OPM_THROW(std::runtime_error, - "The table specified by the PVTW keyword is required" - "to exhibit exactly one row."); - ref_press_ = pvtwTable.getPressureColumn()[0]; - ref_B_ = pvtwTable.getFormationFactorColumn()[0]; - comp_ = pvtwTable.getCompressibilityColumn()[0]; - viscosity_ = pvtwTable.getViscosityColumn()[0]; - visc_comp_ = pvtwTable.getViscosibilityColumn()[0]; - } - - SinglePvtConstCompr(const Opm::PvcdoTable &pvcdoTable) - { - if (pvcdoTable.numRows() != 1) - OPM_THROW(std::runtime_error, - "The table specified by the PVCDO keyword is required" - "to exhibit exactly one row."); - ref_press_ = pvcdoTable.getPressureColumn()[0]; - ref_B_ = pvcdoTable.getFormationFactorColumn()[0]; - comp_ = pvcdoTable.getCompressibilityColumn()[0]; - viscosity_ = pvcdoTable.getViscosityColumn()[0]; - visc_comp_ = pvcdoTable.getViscosibilityColumn()[0]; - } - - SinglePvtConstCompr(double visc) - : ref_press_(0.0), - ref_B_(1.0), - comp_(0.0), - viscosity_(visc), - visc_comp_(0.0) - { - } - - virtual ~SinglePvtConstCompr() - { - } - - virtual void mu(const int n, - const double* p, - const double* /*z*/, - double* output_mu) const - { - if (visc_comp_) { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - double x = -visc_comp_*(p[i] - ref_press_); - output_mu[i] = viscosity_/(1.0 + x + 0.5*x*x); - } - } else { - std::fill(output_mu, output_mu + n, viscosity_); - } - } - - virtual void mu(const int n, - const double* p, - const double* /*r*/, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - if (visc_comp_) { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - double x = -visc_comp_*(p[i] - ref_press_); - double d = (1.0 + x + 0.5*x*x); - output_mu[i] = viscosity_/d; - output_dmudp[i] = (viscosity_/(d*d))*(1+x) * visc_comp_; - } - } else { - std::fill(output_mu, output_mu + n, viscosity_); - std::fill(output_dmudp, output_dmudp + n, 0.0); - } - std::fill(output_dmudr, output_dmudr + n, 0.0); - } - - virtual void mu(const int n, - const double* p, - const double* /*r*/, - const PhasePresence* /*cond*/, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - if (visc_comp_) { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - double x = -visc_comp_*(p[i] - ref_press_); - double d = (1.0 + x + 0.5*x*x); - output_mu[i] = viscosity_/d; - output_dmudp[i] = (viscosity_/(d*d))*(1+x) * visc_comp_; - } - } else { - std::fill(output_mu, output_mu + n, viscosity_); - std::fill(output_dmudp, output_dmudp + n, 0.0); - } - std::fill(output_dmudr, output_dmudr + n, 0.0); - } - - virtual void B(const int n, - const double* p, - const double* /*z*/, - double* output_B) const - { - if (comp_) { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - double x = comp_*(p[i] - ref_press_); - output_B[i] = ref_B_/(1.0 + x + 0.5*x*x); - } - } else { - std::fill(output_B, output_B + n, ref_B_); - } - } - - virtual void dBdp(const int n, - const double* p, - const double* /*z*/, - double* output_B, - double* output_dBdp) const - { - if (comp_) { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - double x = comp_*(p[i] - ref_press_); - double d = (1.0 + x + 0.5*x*x); - output_B[i] = ref_B_/d; - output_dBdp[i] = (-ref_B_/(d*d))*(1 + x) * comp_; - } - } else { - std::fill(output_B, output_B + n, ref_B_); - std::fill(output_dBdp, output_dBdp + n, 0.0); - } - } - - virtual void b(const int n, - const double* p, - const double* /*r*/, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - { - if (comp_) { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - double x = comp_*(p[i] - ref_press_); - double d = (1.0 + x + 0.5*x*x); - - // b = 1/B = d/ref_B_B; - output_b[i] = d/ref_B_; - output_dbdp[i] = (1 + x) * comp_/ref_B_; - } - } else { - - std::fill(output_b, output_b + n, 1/ref_B_); - std::fill(output_dbdp, output_dbdp + n, 0.0); - } - std::fill(output_dbdr, output_dbdr + n, 0.0); - - } - - virtual void b(const int n, - const double* p, - const double* /*r*/, - const PhasePresence* /*cond*/, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - { - if (comp_) { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - double x = comp_*(p[i] - ref_press_); - double d = (1.0 + x + 0.5*x*x); - - // b = 1/B = d/ref_B_B; - output_b[i] = d/ref_B_; - output_dbdp[i] = (1 + x) * comp_/ref_B_; - } - } else { - - std::fill(output_b, output_b + n, 1/ref_B_); - std::fill(output_dbdp, output_dbdp + n, 0.0); - } - std::fill(output_dbdr, output_dbdr + n, 0.0); - - } - - virtual void rsSat(const int n, - const double* /*p*/, - double* output_rsSat, - double* output_drsSatdp) const - - { - std::fill(output_rsSat, output_rsSat + n, 0.0); - std::fill(output_drsSatdp, output_drsSatdp + n, 0.0); - } - - virtual void rvSat(const int n, - const double* /*p*/, - double* output_rvSat, - double* output_drvSatdp) const - - { - std::fill(output_rvSat, output_rvSat + n, 0.0); - std::fill(output_drvSatdp, output_drvSatdp + n, 0.0); - } - - virtual void R(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_R) const - { - std::fill(output_R, output_R + n, 0.0); - } - - virtual void dRdp(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_R, - double* output_dRdp) const - { - std::fill(output_R, output_R + n, 0.0); - std::fill(output_dRdp, output_dRdp + n, 0.0); - } - - private: - double ref_press_; - double ref_B_; - double comp_; - double viscosity_; - double visc_comp_; - }; - -} - - -#endif // OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED - diff --git a/opm/core/props/pvt/SinglePvtDeadSpline.cpp b/opm/core/props/pvt/SinglePvtDeadSpline.cpp deleted file mode 100644 index 9f64197d..00000000 --- a/opm/core/props/pvt/SinglePvtDeadSpline.cpp +++ /dev/null @@ -1,226 +0,0 @@ -/* - Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include -#include - -#include - -#include - -// Extra includes for debug dumping of tables. -// #include -// #include -// #include - -namespace Opm -{ - //------------------------------------------------------------------------ - // Member functions - //------------------------------------------------------------------------- - - SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples) - { - int numRows = pvdoTable.numRows(); - - // Copy data - const std::vector &press = pvdoTable.getPressureColumn(); - const std::vector &B = pvdoTable.getFormationFactorColumn(); - const std::vector &visc = pvdoTable.getViscosityColumn(); - - std::vector B_inv(numRows); - for (int i = 0; i < numRows; ++i) { - B_inv[i] = 1.0 / B[i]; - } - - buildUniformMonotoneTable(press, B_inv, samples, b_); - buildUniformMonotoneTable(press, visc, samples, viscosity_); - } - - SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples) - { - int numRows = pvdgTable.numRows(); - - // Copy data - const std::vector &press = pvdgTable.getPressureColumn(); - const std::vector &B = pvdgTable.getFormationFactorColumn(); - const std::vector &visc = pvdgTable.getViscosityColumn(); - - std::vector B_inv(numRows); - for (int i = 0; i < numRows; ++i) { - B_inv[i] = 1.0 / B[i]; - } - - buildUniformMonotoneTable(press, B_inv, samples, b_); - buildUniformMonotoneTable(press, visc, samples, viscosity_); - } - - // Destructor - SinglePvtDeadSpline::~SinglePvtDeadSpline() - { - } - - - - void SinglePvtDeadSpline::mu(const int n, - const double* p, - const double* /*z*/, - double* output_mu) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_mu[i] = viscosity_(p[i]); - } - } - - void SinglePvtDeadSpline::mu(const int n, - const double* p, - const double* /*r*/, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - // #pragma omp parallel for - - for (int i = 0; i < n; ++i) { - output_mu[i] = viscosity_(p[i]); - output_dmudp[i] = viscosity_.derivative(p[i]); - } - std::fill(output_dmudr, output_dmudr + n, 0.0); - - } - - void SinglePvtDeadSpline::mu(const int n, - const double* p, - const double* /*r*/, - const PhasePresence* /*cond*/, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - // #pragma omp parallel for - - for (int i = 0; i < n; ++i) { - output_mu[i] = viscosity_(p[i]); - output_dmudp[i] = viscosity_.derivative(p[i]); - } - std::fill(output_dmudr, output_dmudr + n, 0.0); - - } - - void SinglePvtDeadSpline::B(const int n, - const double* p, - const double* /*z*/, - double* output_B) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_B[i] = 1.0/b_(p[i]); - } - } - - void SinglePvtDeadSpline::dBdp(const int n, - const double* p, - const double* /*z*/, - double* output_B, - double* output_dBdp) const - { - B(n, p, 0, output_B); -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - double Bg = output_B[i]; - output_dBdp[i] = -Bg*Bg*b_.derivative(p[i]); - } - } - - void SinglePvtDeadSpline::b(const int n, - const double* p, - const double* /*r*/, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_b[i] = b_(p[i]); - output_dbdp[i] = b_.derivative(p[i]); - - } - std::fill(output_dbdr, output_dbdr + n, 0.0); - - } - - void SinglePvtDeadSpline::b(const int n, - const double* p, - const double* /*r*/, - const PhasePresence* /*cond*/, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_b[i] = b_(p[i]); - output_dbdp[i] = b_.derivative(p[i]); - - } - std::fill(output_dbdr, output_dbdr + n, 0.0); - - } - - void SinglePvtDeadSpline::rsSat(const int n, - const double* /*p*/, - double* output_rsSat, - double* output_drsSatdp) const - { - std::fill(output_rsSat, output_rsSat + n, 0.0); - std::fill(output_drsSatdp, output_drsSatdp + n, 0.0); - } - - void SinglePvtDeadSpline::rvSat(const int n, - const double* /*p*/, - double* output_rvSat, - double* output_drvSatdp) const - { - std::fill(output_rvSat, output_rvSat + n, 0.0); - std::fill(output_drvSatdp, output_drvSatdp + n, 0.0); - } - - void SinglePvtDeadSpline::R(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_R) const - { - std::fill(output_R, output_R + n, 0.0); - } - - void SinglePvtDeadSpline::dRdp(const int n, - const double* /*p*/, - const double* /*z*/, - double* output_R, - double* output_dRdp) const - { - std::fill(output_R, output_R + n, 0.0); - std::fill(output_dRdp, output_dRdp + n, 0.0); - } - -} diff --git a/opm/core/props/pvt/SinglePvtInterface.cpp b/opm/core/props/pvt/SinglePvtInterface.cpp deleted file mode 100644 index bbd905e6..00000000 --- a/opm/core/props/pvt/SinglePvtInterface.cpp +++ /dev/null @@ -1,48 +0,0 @@ -/* - Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include - - -namespace Opm -{ - - SinglePvtInterface::SinglePvtInterface() - : num_phases_(MaxNumPhases) - { - for (int i = 0; i < MaxNumPhases; ++i) { - phase_pos_[i] = i; - } - } - - SinglePvtInterface::~SinglePvtInterface() - { - } - - void SinglePvtInterface::setPhaseConfiguration(const int num_phases, const int* phase_pos) - { - num_phases_ = num_phases; - for (int i = 0; i < MaxNumPhases; ++i) { - phase_pos_[i] = phase_pos[i]; - } - } - - -} // namespace Opm diff --git a/opm/core/props/pvt/SinglePvtLiveGas.cpp b/opm/core/props/pvt/SinglePvtLiveGas.cpp deleted file mode 100644 index abc9a641..00000000 --- a/opm/core/props/pvt/SinglePvtLiveGas.cpp +++ /dev/null @@ -1,492 +0,0 @@ -//=========================================================================== -// -// File: MiscibilityLiveGas.cpp -// -// Created: Wed Feb 10 09:21:53 2010 -// -// Author: Bjørn Spjelkavik -// -// Revision: $Id$ -// -//=========================================================================== -/* - Copyright 2010 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include -#include -#include - -#include -#include - -namespace Opm -{ - - using Opm::linearInterpolation; - using Opm::linearInterpolationDerivative; - - - //------------------------------------------------------------------------ - // Member functions - //------------------------------------------------------------------------- - SinglePvtLiveGas::SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable) - { - // GAS, PVTG - saturated_gas_table_.resize(4); - saturated_gas_table_[0] = pvtgTable.getOuterTable()->getPressureColumn(); - saturated_gas_table_[1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); - saturated_gas_table_[2] = pvtgTable.getOuterTable()->getGasViscosityColumn(); - saturated_gas_table_[3] = pvtgTable.getOuterTable()->getOilSolubilityColumn(); - - int sz = pvtgTable.getOuterTable()->numRows(); - undersat_gas_tables_.resize(sz); - for (int i=0; i 1/Bg - for (int i=0; i saturated_gas_table_[0][ltp]) { - return linearInterpolation(undersat_gas_tables_[ltp][0], - undersat_gas_tables_[ltp][item], - maxR); - } - - // Interpolate between table sections - double w = (press - saturated_gas_table_[0][is]) / - (saturated_gas_table_[0][is+1] - - saturated_gas_table_[0][is]); - if (undersat_gas_tables_[is][0].size() < 2) { - double val = saturated_gas_table_[item][is] + - w*(saturated_gas_table_[item][is+1] - - saturated_gas_table_[item][is]); - return val; - } - double val1 = - linearInterpolation(undersat_gas_tables_[is][0], - undersat_gas_tables_[is][item], - maxR); - double val2 = - linearInterpolation(undersat_gas_tables_[is+1][0], - undersat_gas_tables_[is+1][item], - maxR); - double val = val1 + w*(val2 - val1); - return val; - } - } - } - - double SinglePvtLiveGas::miscible_gas(const double press, - const double r, - const PhasePresence& cond, - const int item, - const int deriv) const -{ - const bool isSat = cond.hasFreeOil(); - - // Derivative w.r.t p - if (deriv == 1) { - if (isSat) { // Saturated case - return linearInterpolationDerivative(saturated_gas_table_[0], - saturated_gas_table_[item], - press); - } else { // Undersaturated case - int is = tableIndex(saturated_gas_table_[0], press); - if (undersat_gas_tables_[is][0].size() < 2) { - double val = (saturated_gas_table_[item][is+1] - - saturated_gas_table_[item][is]) / - (saturated_gas_table_[0][is+1] - - saturated_gas_table_[0][is]); - return val; - } - double val1 = - linearInterpolation(undersat_gas_tables_[is][0], - undersat_gas_tables_[is][item], - r); - double val2 = - linearInterpolation(undersat_gas_tables_[is+1][0], - undersat_gas_tables_[is+1][item], - r); - double val = (val2 - val1)/ - (saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]); - return val; - } - } else if (deriv == 2){ - if (isSat) { - return 0; - } else { - int is = tableIndex(saturated_gas_table_[0], press); - double w = (press - saturated_gas_table_[0][is]) / - (saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]); - assert(undersat_gas_tables_[is][0].size() >= 2); - assert(undersat_gas_tables_[is+1][0].size() >= 2); - double val1 = - linearInterpolationDerivative(undersat_gas_tables_[is][0], - undersat_gas_tables_[is][item], - r); - double val2 = - linearInterpolationDerivative(undersat_gas_tables_[is+1][0], - undersat_gas_tables_[is+1][item], - r); - - double val = val1 + w * (val2 - val1); - return val; - - } - } else { - if (isSat) { // Saturated case - return linearInterpolation(saturated_gas_table_[0], - saturated_gas_table_[item], - press); - } else { // Undersaturated case - int is = tableIndex(saturated_gas_table_[0], press); - // Extrapolate from first table section - if (is == 0 && press < saturated_gas_table_[0][0]) { - return linearInterpolation(undersat_gas_tables_[0][0], - undersat_gas_tables_[0][item], - r); - } - - // Extrapolate from last table section - //int ltp = saturated_gas_table_[0].size() - 1; - //if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) { - // return linearInterpolation(undersat_gas_tables_[ltp][0], - // undersat_gas_tables_[ltp][item], - // r); - //} - - // Interpolate between table sections - double w = (press - saturated_gas_table_[0][is]) / - (saturated_gas_table_[0][is+1] - - saturated_gas_table_[0][is]); - if (undersat_gas_tables_[is][0].size() < 2) { - double val = saturated_gas_table_[item][is] + - w*(saturated_gas_table_[item][is+1] - - saturated_gas_table_[item][is]); - return val; - } - double val1 = - linearInterpolation(undersat_gas_tables_[is][0], - undersat_gas_tables_[is][item], - r); - double val2 = - linearInterpolation(undersat_gas_tables_[is+1][0], - undersat_gas_tables_[is+1][item], - r); - double val = val1 + w*(val2 - val1); - return val; - } - } - } - - - - -} // namespace Opm diff --git a/opm/core/props/pvt/SinglePvtLiveOil.cpp b/opm/core/props/pvt/SinglePvtLiveOil.cpp deleted file mode 100644 index b7e4b113..00000000 --- a/opm/core/props/pvt/SinglePvtLiveOil.cpp +++ /dev/null @@ -1,551 +0,0 @@ -/* - Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include -#include -#include - -#include - -namespace Opm -{ - - using Opm::linearInterpolation; - using Opm::linearInterpolationDerivative; - using Opm::tableIndex; - - //------------------------------------------------------------------------ - // Member functions - //------------------------------------------------------------------------- - - SinglePvtLiveOil::SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable) - { - const auto saturatedPvto = pvtoTable.getOuterTable(); - - // OIL, PVTO - saturated_oil_table_.resize(4); - const int sz = saturatedPvto->numRows(); - for (int k=0; k<4; ++k) { - saturated_oil_table_[k].resize(sz); - } - for (int i=0; igetPressureColumn()[i]; // p - saturated_oil_table_[1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo - saturated_oil_table_[2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o - saturated_oil_table_[3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs - } - - undersat_oil_tables_.resize(sz); - for (int i=0; inumRows(); - undersat_oil_tables_[i][0].resize(tsize); - undersat_oil_tables_[i][1].resize(tsize); - undersat_oil_tables_[i][2].resize(tsize); - for (int j=0; jgetPressureColumn()[j]; // p - undersat_oil_tables_[i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j]; // 1/Bo - undersat_oil_tables_[i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j]; // mu_o - } - } - - // Complete undersaturated tables by extrapolating from existing data - // as is done in Eclipse and Mrst - int iNext = -1; - for (int i=0; i 1) { - continue; - } - // Look ahead for next record containing undersaturated data - if (iNext < i) { - iNext = i+1; - while (iNext > >::size_type sz_t; - for (sz_t j=1; j= 2); - assert(undersat_oil_tables_[is+1][0].size() >= 2); - double val1 = - linearInterpolationDerivative(undersat_oil_tables_[is][0], - undersat_oil_tables_[is][item], - press); - double val2 = - linearInterpolationDerivative(undersat_oil_tables_[is+1][0], - undersat_oil_tables_[is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - } else { - if (Rval < maxR ) { // Saturated case - return linearInterpolation(saturated_oil_table_[0], - saturated_oil_table_[item], - press); - } else { // Undersaturated case - // Interpolate between table sections - int is = tableIndex(saturated_oil_table_[3], maxR); - double w = (maxR - saturated_oil_table_[3][is]) / - (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]); - assert(undersat_oil_tables_[is][0].size() >= 2); - assert(undersat_oil_tables_[is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[is][0], - undersat_oil_tables_[is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[is+1][0], - undersat_oil_tables_[is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - } - } - - double SinglePvtLiveOil::miscible_oil(const double press, - const double r, - const int item, - const int deriv) const - { - int section; - double Rval = linearInterpolation(saturated_oil_table_[0], - saturated_oil_table_[3], - press, section); - // derivative with respect to frist component (pressure) - if (deriv == 1) { - if (Rval <= r ) { // Saturated case - return linearInterpolationDerivative(saturated_oil_table_[0], - saturated_oil_table_[item], - press); - } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[3], r); - double w = (r - saturated_oil_table_[3][is]) / - (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]); - assert(undersat_oil_tables_[is][0].size() >= 2); - assert(undersat_oil_tables_[is+1][0].size() >= 2); - double val1 = - linearInterpolationDerivative(undersat_oil_tables_[is][0], - undersat_oil_tables_[is][item], - press); - double val2 = - linearInterpolationDerivative(undersat_oil_tables_[is+1][0], - undersat_oil_tables_[is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - // derivative with respect to second component (r) - } else if (deriv == 2) { - if (Rval <= r ) { // Saturated case - return 0; - } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[3], r); - assert(undersat_oil_tables_[is][0].size() >= 2); - assert(undersat_oil_tables_[is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[is][0], - undersat_oil_tables_[is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[is+1][0], - undersat_oil_tables_[is+1][item], - press); - - double val = (val2 - val1)/(saturated_oil_table_[3][is+1]-saturated_oil_table_[3][is]); - return val; - } - - - } else { - if (Rval <= r ) { // Saturated case - return linearInterpolation(saturated_oil_table_[0], - saturated_oil_table_[item], - press); - } else { // Undersaturated case - // Interpolate between table sections - int is = tableIndex(saturated_oil_table_[3], r); - double w = (r - saturated_oil_table_[3][is]) / - (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]); - assert(undersat_oil_tables_[is][0].size() >= 2); - assert(undersat_oil_tables_[is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[is][0], - undersat_oil_tables_[is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[is+1][0], - undersat_oil_tables_[is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - } - } - - double SinglePvtLiveOil::miscible_oil(const double press, - const double r, - const PhasePresence& cond, - const int item, - const int deriv) const - { - const bool isSat = cond.hasFreeGas(); - - // derivative with respect to frist component (pressure) - if (deriv == 1) { - if (isSat) { // Saturated case - return linearInterpolationDerivative(saturated_oil_table_[0], - saturated_oil_table_[item], - press); - } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[3], r); - double w = (r - saturated_oil_table_[3][is]) / - (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]); - assert(undersat_oil_tables_[is][0].size() >= 2); - assert(undersat_oil_tables_[is+1][0].size() >= 2); - double val1 = - linearInterpolationDerivative(undersat_oil_tables_[is][0], - undersat_oil_tables_[is][item], - press); - double val2 = - linearInterpolationDerivative(undersat_oil_tables_[is+1][0], - undersat_oil_tables_[is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - // derivative with respect to second component (r) - } else if (deriv == 2) { - if (isSat) { // Saturated case - return 0; - } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[3], r); - assert(undersat_oil_tables_[is][0].size() >= 2); - assert(undersat_oil_tables_[is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[is][0], - undersat_oil_tables_[is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[is+1][0], - undersat_oil_tables_[is+1][item], - press); - - double val = (val2 - val1)/(saturated_oil_table_[3][is+1]-saturated_oil_table_[3][is]); - return val; - } - - - } else { - if (isSat) { // Saturated case - return linearInterpolation(saturated_oil_table_[0], - saturated_oil_table_[item], - press); - } else { // Undersaturated case - // Interpolate between table sections - int is = tableIndex(saturated_oil_table_[3], r); - double w = (r - saturated_oil_table_[3][is]) / - (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]); - assert(undersat_oil_tables_[is][0].size() >= 2); - assert(undersat_oil_tables_[is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[is][0], - undersat_oil_tables_[is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[is+1][0], - undersat_oil_tables_[is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - } - } - -} // namespace Opm diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 1cf88ee4..e5f8c8ec 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -279,7 +279,7 @@ namespace Opm if (rec[i].live_oil_table_index > 0) { if (deck->hasKeyword("RSVD") && size_t(rec[i].live_oil_table_index) <= deck->getKeyword("RSVD")->size()) { - Opm::SingleRecordTable rsvd(deck->getKeyword("RSVD"),std::vector{"vd", "rs"},rec[i].live_oil_table_index-1); + Opm::SingleRecordTable rsvd(deck->getKeyword("RSVD"),std::vector{"vd", "rs"},rec[i].live_oil_table_index-1); std::vector vd(rsvd.getColumn("vd")); std::vector rs(rsvd.getColumn("rs")); rs_func_.push_back(std::make_shared(props, cell, vd, rs)); diff --git a/tests/liveoil.DATA b/tests/liveoil.DATA index 48fcafe6..b3c3e271 100644 --- a/tests/liveoil.DATA +++ b/tests/liveoil.DATA @@ -9,6 +9,19 @@ WATER FIELD +-- tests for the PVT functions need a grid because the OPM-API for the +-- PVT functions assumes the presence of compressed grid cells, +-- i.e. we need to be able to map from compressed to logical cartesian +-- cell indices to find out the PVT region of a cell +DIMENS +3 3 3 / +DXV +1 2 3 / +DYV +1 2 3 / +DZV +1 2 3 / + PVTO -- Rs Pbub Bo Vo .0 14.7 1.0000 1.20 / diff --git a/tests/test_blackoilfluid.cpp b/tests/test_blackoilfluid.cpp index 5e2107f4..8dcf3417 100644 --- a/tests/test_blackoilfluid.cpp +++ b/tests/test_blackoilfluid.cpp @@ -1,15 +1,17 @@ #include #include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include #include #include +#include #include #include +#include #include #include @@ -36,63 +38,78 @@ using namespace Opm; using namespace std; -std::vector > getProps(Opm::DeckConstPtr deck,PhaseUsage phase_usage_){ - +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_; + std::vector > props_; // Set the properties. props_.resize(phase_usage_.num_phases); // Water PVT if (phase_usage_.phase_used[Aqua]) { if (deck->hasKeyword("PVTW")) { - Opm::PvtwTable pvtwTable(deck->getKeyword("PVTW")); - props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable)); + std::shared_ptr pvtw(new PvtConstCompr); + pvtw->initFromWater(deck->getKeyword("PVTW"), pvtTableIdx); + + props_[phase_usage_.phase_pos[Aqua]] = pvtw; } else { // Eclipse 100 default. - props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise)); + props_[phase_usage_.phase_pos[Aqua]].reset(new PvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise)); } } // Oil PVT if (phase_usage_.phase_used[Liquid]) { if (deck->hasKeyword("PVDO")) { - Opm::PvdoTable pvdoTable(deck->getKeyword("PVDO")); + Opm::DeckKeywordConstPtr pvdoKeyword(deck->getKeyword("PVDO")); if (samples > 0) { - props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples)); + std::shared_ptr splinePvt(new PvtDeadSpline); + splinePvt->initFromOil(pvdoKeyword, pvtTableIdx, samples); + props_[phase_usage_.phase_pos[Liquid]] = splinePvt; } else { - props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable)); + std::shared_ptr deadPvt(new PvtDead); + deadPvt->initFromOil(pvdoKeyword, pvtTableIdx); + props_[phase_usage_.phase_pos[Liquid]] = deadPvt; } } else if (deck->hasKeyword("PVTO")) { - Opm::PvtoTable pvtoTable(deck->getKeyword("PVTO"), /*tableIdx=*/0); - - props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable)); + props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"), pvtTableIdx)); } else if (deck->hasKeyword("PVCDO")) { - Opm::PvcdoTable pvcdoTable(deck->getKeyword("PVCDO"), /*tableIdx=*/0); + std::shared_ptr pvcdo(new PvtConstCompr); + pvcdo->initFromOil(deck->getKeyword("PVCDO"), pvtTableIdx); - props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable)); + props_[phase_usage_.phase_pos[Liquid]] = pvcdo; } else { - OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n"); + OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n"); } } // Gas PVT if (phase_usage_.phase_used[Vapour]) { if (deck->hasKeyword("PVDG")) { - Opm::PvdgTable pvdgTable(deck->getKeyword("PVDG")); - + Opm::DeckKeywordConstPtr pvdgKeyword(deck->getKeyword("PVDG")); if (samples > 0) { - props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples)); + std::shared_ptr splinePvt(new PvtDeadSpline); + splinePvt->initFromGas(pvdgKeyword, pvtTableIdx, samples); + props_[phase_usage_.phase_pos[Vapour]] = splinePvt; } else { - props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable)); + std::shared_ptr deadPvt(new PvtDead); + deadPvt->initFromGas(pvdgKeyword, pvtTableIdx); + props_[phase_usage_.phase_pos[Vapour]] = deadPvt; } } else if (deck->hasKeyword("PVTG")) { - Opm::PvtgTable pvtgTable(deck->getKeyword("PVTG"), /*tableIdx=*/0); - - props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable)); + props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG"), pvtTableIdx)); } else { OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n"); } @@ -102,7 +119,7 @@ std::vector > getProps(Opm::DeckConstPtr dec } void testmu(const double reltol, int n, int np, std::vector p, std::vector r,std::vector z, - std::vector > props_, std::vector condition) + std::vector > props_, std::vector condition) { std::vector mu(n); std::vector dmudp(n); @@ -138,7 +155,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, - std::vector > props_, std::vector condition) + std::vector > props_, std::vector condition) { // test b std::vector b(n); @@ -180,7 +197,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, std::vector p, std::vector > props_){ // test bublepoint pressure std::vector rs(n); std::vector drsdp(n); @@ -202,7 +219,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, std::vector p, std::vector > props_){ // test rv saturated std::vector rv(n); std::vector drvdp(n); @@ -234,9 +251,11 @@ BOOST_AUTO_TEST_CASE(test_liveoil) Opm::ParserPtr parser(new Opm::Parser()); Opm::DeckConstPtr deck(parser->parseFile(filename)); + + // setup pvt interface PhaseUsage phase_usage_ = phaseUsageFromDeck(deck); - std::vector > props_ = getProps(deck,phase_usage_); + std::vector > props_ = getProps(deck,phase_usage_); // setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference @@ -311,7 +330,7 @@ BOOST_AUTO_TEST_CASE(test_wetgas) // setup pvt interface PhaseUsage phase_usage_ = phaseUsageFromDeck(deck); - std::vector > props_ = getProps(deck,phase_usage_); + std::vector > props_ = getProps(deck,phase_usage_); // setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference diff --git a/tests/wetgas.DATA b/tests/wetgas.DATA index 3c094dc1..1a98a975 100644 --- a/tests/wetgas.DATA +++ b/tests/wetgas.DATA @@ -9,6 +9,19 @@ GAS FIELD +-- tests for the PVT functions need a grid because the OPM-API for the +-- PVT functions assumes the presence of compressed grid cells, +-- i.e. we need to be able to map from compressed to logical cartesian +-- cell indices to find out the PVT region of a cell +DIMENS +3 3 3 / +DXV +1 2 3 / +DYV +1 2 3 / +DZV +1 2 3 / + -- PVT PROPERTIES OF DRY GAS (NO VAPOURISED OIL) -- FROM SPE3 Blackoil Kleppe --