diff --git a/opm/autodiff/AutoDiffHelpers.hpp b/opm/autodiff/AutoDiffHelpers.hpp index 00672a348..bbd631c0e 100644 --- a/opm/autodiff/AutoDiffHelpers.hpp +++ b/opm/autodiff/AutoDiffHelpers.hpp @@ -66,9 +66,7 @@ struct HelperOps const int nc = numCells(grid); const int nf = numFaces(grid); // Define some neighbourhood-derived helper arrays. - typedef Eigen::Array OneColBool; typedef Eigen::Array TwoColInt; - typedef Eigen::Array TwoColBool; TwoColInt nbi; extractInternalFaces(grid, internal_faces, nbi); int num_internal=internal_faces.size(); diff --git a/opm/autodiff/BlackoilPropsAd.cpp b/opm/autodiff/BlackoilPropsAd.cpp index 3371c5ae0..22288afff 100644 --- a/opm/autodiff/BlackoilPropsAd.cpp +++ b/opm/autodiff/BlackoilPropsAd.cpp @@ -91,8 +91,10 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of 3 density values. - const double* BlackoilPropsAd::surfaceDensity() const + const double* BlackoilPropsAd::surfaceDensity(int regionIdx) const { + // this class only supports a single PVT region for now... + assert(regionIdx == 0); return props_.surfaceDensity(); } diff --git a/opm/autodiff/BlackoilPropsAd.hpp b/opm/autodiff/BlackoilPropsAd.hpp index 5a8f35e98..b972ab76c 100644 --- a/opm/autodiff/BlackoilPropsAd.hpp +++ b/opm/autodiff/BlackoilPropsAd.hpp @@ -94,7 +94,7 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of 3 density values. - const double* surfaceDensity() const; + const double* surfaceDensity(int regionIdx = 0) const; // ------ Viscosity ------ diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.cpp b/opm/autodiff/BlackoilPropsAdFromDeck.cpp index 554ef8305..698f3dc39 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.cpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.cpp @@ -23,12 +23,12 @@ #include #include #include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include #include #include @@ -80,79 +80,104 @@ namespace Opm int dimension, const bool init_rock) { + // retrieve the cell specific PVT table index from the deck + // and using the grid... + extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell); + if (init_rock){ rock_.init(deck, number_of_cells, global_cell, cart_dims); } - const int samples = 0; - const int region_number = 0; - phase_usage_ = phaseUsageFromDeck(deck); // Surface densities. Accounting for different orders in eclipse and our code. - if (deck->hasKeyword("DENSITY")) { - const auto keyword = deck->getKeyword("DENSITY"); - const auto record = keyword->getRecord(region_number); - enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 }; + 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_[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]] = record->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]] = record->getItem("GAS")->getSIDouble(0); + densities_[regionIdx][phase_usage_.phase_pos[Vapour]] + = densityKeyword->getRecord(regionIdx)->getItem("GAS")->getSIDouble(0); } - if (phase_usage_.phase_used[Liquid]) { - densities_[phase_usage_.phase_pos[Liquid]] = record->getItem("OIL")->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, + number_of_cells, + global_cell); + + const int numSamples = 0; + + // 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"), region_number); - 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)); - } - } + // 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")); + 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, 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); + 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)); + } else if (deck->hasKeyword("PVTO")) { + props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"))); } else if (deck->hasKeyword("PVCDO")) { - Opm::PvcdoTable pvdcoTable(deck->getKeyword("PVCDO"), region_number); - props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvdcoTable)); + std::shared_ptr pvcdo(new PvtConstCompr); + pvcdo->initFromOil(deck->getKeyword("PVCDO")); + + props_[phase_usage_.phase_pos[Liquid]] = pvcdo; } else { - OPM_THROW(std::runtime_error, "Input is missing PVDO, PVTO or PVCDO\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::PvdoTable 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, 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); + + 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"))); } else { OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n"); } @@ -221,9 +246,10 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of 3 density values. - const double* BlackoilPropsAdFromDeck::surfaceDensity() const + const double* BlackoilPropsAdFromDeck::surfaceDensity(const int cellIdx) const { - return densities_; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + return &densities_[pvtRegionIdx][0]; } @@ -246,7 +272,7 @@ namespace Opm V dmudr(n); const double* rs = 0; - props_[phase_usage_.phase_pos[Water]]->mu(n, pw.data(), rs, + props_[phase_usage_.phase_pos[Water]]->mu(n, &pvtTableIdx_[0], pw.data(), rs, mu.data(), dmudp.data(), dmudr.data()); return mu; } @@ -271,7 +297,7 @@ namespace Opm V dmudp(n); V dmudr(n); - props_[phase_usage_.phase_pos[Oil]]->mu(n, po.data(), rs.data(), &cond[0], + props_[phase_usage_.phase_pos[Oil]]->mu(n, &pvtTableIdx_[0], po.data(), rs.data(), &cond[0], mu.data(), dmudp.data(), dmudr.data()); return mu; } @@ -293,7 +319,7 @@ namespace Opm V dmudr(n); const double* rs = 0; - props_[phase_usage_.phase_pos[Gas]]->mu(n, pg.data(), rs, + props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.data(), rs, mu.data(), dmudp.data(), dmudr.data()); return mu; } @@ -316,7 +342,7 @@ namespace Opm V dmudp(n); V dmudr(n); - props_[phase_usage_.phase_pos[Gas]]->mu(n, pg.data(), rv.data(),&cond[0], + props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.data(), rv.data(),&cond[0], mu.data(), dmudp.data(), dmudr.data()); return mu; } @@ -338,7 +364,7 @@ namespace Opm V dmudr(n); const double* rs = 0; - props_[phase_usage_.phase_pos[Water]]->mu(n, pw.value().data(), rs, + props_[phase_usage_.phase_pos[Water]]->mu(n, &pvtTableIdx_[0], pw.value().data(), rs, mu.data(), dmudp.data(), dmudr.data()); ADB::M dmudp_diag = spdiag(dmudp); const int num_blocks = pw.numBlocks(); @@ -369,7 +395,7 @@ namespace Opm V dmudp(n); V dmudr(n); - props_[phase_usage_.phase_pos[Oil]]->mu(n, po.value().data(), rs.value().data(), + props_[phase_usage_.phase_pos[Oil]]->mu(n, &pvtTableIdx_[0], po.value().data(), rs.value().data(), &cond[0], mu.data(), dmudp.data(), dmudr.data()); ADB::M dmudp_diag = spdiag(dmudp); @@ -399,7 +425,7 @@ namespace Opm V dmudr(n); const double* rv = 0; - props_[phase_usage_.phase_pos[Gas]]->mu(n, pg.value().data(), rv, + props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.value().data(), rv, mu.data(), dmudp.data(), dmudr.data()); ADB::M dmudp_diag = spdiag(dmudp); @@ -431,7 +457,7 @@ namespace Opm V dmudp(n); V dmudr(n); - props_[phase_usage_.phase_pos[Gas]]->mu(n, pg.value().data(), rv.value().data(),&cond[0], + props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.value().data(), rv.value().data(),&cond[0], mu.data(), dmudp.data(), dmudr.data()); ADB::M dmudp_diag = spdiag(dmudp); @@ -480,7 +506,7 @@ namespace Opm V dbdr(n); const double* rs = 0; - props_[phase_usage_.phase_pos[Water]]->b(n, pw.data(), rs, + props_[phase_usage_.phase_pos[Water]]->b(n, &pvtTableIdx_[0], pw.data(), rs, b.data(), dbdp.data(), dbdr.data()); return b; @@ -507,7 +533,7 @@ namespace Opm V dbdp(n); V dbdr(n); - props_[phase_usage_.phase_pos[Oil]]->b(n, po.data(), rs.data(), &cond[0], + props_[phase_usage_.phase_pos[Oil]]->b(n, &pvtTableIdx_[0], po.data(), rs.data(), &cond[0], b.data(), dbdp.data(), dbdr.data()); return b; @@ -531,7 +557,7 @@ namespace Opm V dbdr(n); const double* rs = 0; - props_[phase_usage_.phase_pos[Gas]]->b(n, pg.data(), rs, + props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.data(), rs, b.data(), dbdp.data(), dbdr.data()); return b; @@ -558,7 +584,7 @@ namespace Opm V dbdp(n); V dbdr(n); - props_[phase_usage_.phase_pos[Gas]]->b(n, pg.data(), rv.data(), &cond[0], + props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.data(), rv.data(), &cond[0], b.data(), dbdp.data(), dbdr.data()); return b; @@ -582,7 +608,7 @@ namespace Opm V dbdr(n); const double* rs = 0; - props_[phase_usage_.phase_pos[Water]]->b(n, pw.value().data(), rs, + props_[phase_usage_.phase_pos[Water]]->b(n, &pvtTableIdx_[0], pw.value().data(), rs, b.data(), dbdp.data(), dbdr.data()); ADB::M dbdp_diag = spdiag(dbdp); @@ -615,7 +641,7 @@ namespace Opm V dbdp(n); V dbdr(n); - props_[phase_usage_.phase_pos[Oil]]->b(n, po.value().data(), rs.value().data(), + props_[phase_usage_.phase_pos[Oil]]->b(n, &pvtTableIdx_[0], po.value().data(), rs.value().data(), &cond[0], b.data(), dbdp.data(), dbdr.data()); ADB::M dbdp_diag = spdiag(dbdp); @@ -646,7 +672,7 @@ namespace Opm V dbdr(n); const double* rv = 0; - props_[phase_usage_.phase_pos[Gas]]->b(n, pg.value().data(), rv, + props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.value().data(), rv, b.data(), dbdp.data(), dbdr.data()); ADB::M dbdp_diag = spdiag(dbdp); @@ -679,7 +705,7 @@ namespace Opm V dbdp(n); V dbdr(n); - props_[phase_usage_.phase_pos[Gas]]->b(n, pg.value().data(), rv.value().data(), &cond[0], + props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.value().data(), rv.value().data(), &cond[0], b.data(), dbdp.data(), dbdr.data()); ADB::M dbdp_diag = spdiag(dbdp); @@ -710,7 +736,7 @@ namespace Opm assert(po.size() == n); V rbub(n); V drbubdp(n); - props_[Oil]->rsSat(n, po.data(), rbub.data(), drbubdp.data()); + props_[Oil]->rsSat(n, &pvtTableIdx_[0], po.data(), rbub.data(), drbubdp.data()); return rbub; } @@ -728,7 +754,7 @@ namespace Opm assert(po.size() == n); V rbub(n); V drbubdp(n); - props_[Oil]->rsSat(n, po.value().data(), rbub.data(), drbubdp.data()); + props_[Oil]->rsSat(n, &pvtTableIdx_[0], po.value().data(), rbub.data(), drbubdp.data()); ADB::M drbubdp_diag = spdiag(drbubdp); const int num_blocks = po.numBlocks(); std::vector jacs(num_blocks); @@ -754,7 +780,7 @@ namespace Opm assert(po.size() == n); V rv(n); V drvdp(n); - props_[Gas]->rvSat(n, po.data(), rv.data(), drvdp.data()); + props_[Gas]->rvSat(n, &pvtTableIdx_[0], po.data(), rv.data(), drvdp.data()); return rv; } @@ -772,7 +798,7 @@ namespace Opm assert(po.size() == n); V rv(n); V drvdp(n); - props_[Gas]->rvSat(n, po.value().data(), rv.data(), drvdp.data()); + props_[Gas]->rvSat(n, &pvtTableIdx_[0], po.value().data(), rv.data(), drvdp.data()); ADB::M drvdp_diag = spdiag(drvdp); const int num_blocks = po.numBlocks(); std::vector jacs(num_blocks); diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.hpp b/opm/autodiff/BlackoilPropsAdFromDeck.hpp index 7478c3887..a491a7e9e 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.hpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.hpp @@ -23,12 +23,14 @@ #include #include +#include #include #include #include #include +#include #ifdef HAVE_DUNE_CORNERPOINT #include "disable_warning_pragmas.h" @@ -38,8 +40,7 @@ namespace Opm { - - class SinglePvtInterface; + class PvtInterface; /// This class implements the AD-adapted fluid interface for /// three-phase black-oil. It requires an input deck from which it @@ -76,6 +77,11 @@ namespace Opm /// \return N, the number of cells. int numCells() const; + /// Return an array containing the PVT table index for each + /// grid cell + virtual const int* cellPvtRegionIndex() const + { return &cellPvtRegionIdx_[0]; } + /// \return Array of N porosity values. const double* porosity() const; @@ -109,7 +115,7 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of 3 density values. - const double* surfaceDensity() const; + const double* surfaceDensity(const int cellIdx = 0) const; // ------ Viscosity ------ @@ -347,12 +353,27 @@ namespace Opm RockFromDeck rock_; std::unique_ptr satprops_; + PhaseUsage phase_usage_; - std::vector > props_; - double densities_[BlackoilPhases::MaxNumPhases]; + + // The PVT region which is to be used for each cell + std::vector cellPvtRegionIdx_; + + // The PVT properties. One object per PVT region and per + // active fluid phase. + std::vector > props_; + + // The index of the PVT table which ought to be used for each + // cell. Eclipse does not seem to allow specifying fluid-phase + // specific table indices, so for the sake of simplicity, we + // don't do this either. (if it turns out that Eclipes does in + // fact support it or if it by some miracle gains this feature + // in the future, this attribute needs to become a vector of + // vectors of ints.) + std::vector pvtTableIdx_; + + std::vector > densities_; }; - - } // namespace Opm #endif // OPM_BLACKOILPROPSADFROMDECK_HEADER_INCLUDED diff --git a/opm/autodiff/BlackoilPropsAdInterface.hpp b/opm/autodiff/BlackoilPropsAdInterface.hpp index d86560543..bde9f2b7f 100644 --- a/opm/autodiff/BlackoilPropsAdInterface.hpp +++ b/opm/autodiff/BlackoilPropsAdInterface.hpp @@ -84,7 +84,7 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of 3 density values. - virtual const double* surfaceDensity() const = 0; + virtual const double* surfaceDensity(int regionIdx = 0) const = 0; // ------ Viscosity ------