diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 606e04bab..24ac6438c 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -238,7 +238,7 @@ namespace Opm const int cell = wells_->well_cells[j]; const double cell_depth = grid_.cell_centroids[dim * cell + dim - 1]; props_.matrix(1, &state.pressure()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0); - props_.density(1, &A[0], &rho[0]); + props_.density(1, &A[0], &cell, &rho[0]); for (int phase = 0; phase < np; ++phase) { const double s_phase = state.saturation()[np*cell + phase]; wellperf_wdp_[j] += s_phase*rho[phase]*grav*(cell_depth - ref_depth); @@ -380,7 +380,7 @@ namespace Opm // Gravity contribution, gravcontrib = rho*(face_z - cell_z) [per phase]. if (grav != 0.0) { const double depth_diff = face_depth - grid_.cell_centroids[c[j]*dim + dim - 1]; - props_.density(1, &cell_A_[np*np*c[j]], &gravcontrib[j][0]); + props_.density(1, &cell_A_[np*np*c[j]], &c[j], &gravcontrib[j][0]); for (int p = 0; p < np; ++p) { gravcontrib[j][p] *= depth_diff*grav; } diff --git a/opm/core/props/BlackoilPropertiesBasic.cpp b/opm/core/props/BlackoilPropertiesBasic.cpp index 8700af117..ab9b238bc 100644 --- a/opm/core/props/BlackoilPropertiesBasic.cpp +++ b/opm/core/props/BlackoilPropertiesBasic.cpp @@ -157,9 +157,11 @@ namespace Opm /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices /// are assumed to be in Fortran order, and are typically the result /// of a call to the method matrix(). + /// \param[in] cells The index of the grid cell of each data point. /// \param[out] rho Array of nP density values, array must be valid before calling. void BlackoilPropertiesBasic::density(const int n, const double* A, + const int* /*cells*/, double* rho) const { const int np = numPhases(); @@ -177,7 +179,7 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of P density values. - const double* BlackoilPropertiesBasic::surfaceDensity() const + const double* BlackoilPropertiesBasic::surfaceDensity(int /*cellIdx*/) const { return pvt_.surfaceDensities(); } diff --git a/opm/core/props/BlackoilPropertiesBasic.hpp b/opm/core/props/BlackoilPropertiesBasic.hpp index ad84b1a32..b85991eda 100644 --- a/opm/core/props/BlackoilPropertiesBasic.hpp +++ b/opm/core/props/BlackoilPropertiesBasic.hpp @@ -59,6 +59,11 @@ namespace Opm /// \return N, the number of cells. virtual int numCells() const; + /// Return an array containing the PVT table index for each + /// grid cell + virtual const int* cellPvtRegionIndex() const + { return NULL; } + /// \return Array of N porosity values. virtual const double* porosity() const; @@ -114,14 +119,17 @@ namespace Opm /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices /// are assumed to be in Fortran order, and are typically the result /// of a call to the method matrix(). + /// \param[in] cells The index of the grid cell of each data point. /// \param[out] rho Array of nP density values, array must be valid before calling. virtual void density(const int n, const double* A, + const int* cells, double* rho) const; /// Densities of stock components at surface conditions. + /// \param[in] cellIdx The index of the cell for which the surface density ought to be calculated /// \return Array of P density values. - virtual const double* surfaceDensity() const; + virtual const double* surfaceDensity(int cellIdx = 0) const; /// \param[in] n Number of data points. /// \param[in] s Array of nP saturation values. diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index 5a781031f..71246f9c2 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -96,14 +96,19 @@ namespace Opm void BlackoilPropertiesFromDeck::viscosity(const int n, const double* p, const double* z, - const int* /*cells*/, + const int* cells, double* mu, double* dmudp) const { if (dmudp) { OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::viscosity() -- derivatives of viscosity not yet implemented."); } else { - pvt_.mu(n, p, z, mu); + const int *cellPvtTableIdx = cellPvtRegionIndex(); + std::vector pvtTableIdx(n); + for (int i = 0; i < n; ++ i) + pvtTableIdx[i] = cellPvtTableIdx[cells[i]]; + + pvt_.mu(n, &pvtTableIdx[0], p, z, mu); } } @@ -120,21 +125,27 @@ namespace Opm void BlackoilPropertiesFromDeck::matrix(const int n, const double* p, const double* z, - const int* /*cells*/, + const int* cells, double* A, double* dAdp) const { const int np = numPhases(); + + const int *cellPvtTableIdx = cellPvtRegionIndex(); + std::vector pvtTableIdx(n); + for (int i = 0; i < n; ++ i) + pvtTableIdx[i] = cellPvtTableIdx[cells[i]]; + B_.resize(n*np); R_.resize(n*np); if (dAdp) { dB_.resize(n*np); dR_.resize(n*np); - pvt_.dBdp(n, p, z, &B_[0], &dB_[0]); - pvt_.dRdp(n, p, z, &R_[0], &dR_[0]); + pvt_.dBdp(n, &pvtTableIdx[0], p, z, &B_[0], &dB_[0]); + pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]); } else { - pvt_.B(n, p, z, &B_[0]); - pvt_.R(n, p, z, &R_[0]); + pvt_.B(n, &pvtTableIdx[0], p, z, &B_[0]); + pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]); } const int* phase_pos = pvt_.phasePosition(); bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] && @@ -207,15 +218,19 @@ namespace Opm /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices /// are assumed to be in Fortran order, and are typically the result /// of a call to the method matrix(). + /// \param[in] cells The index of the grid cell of each data point. /// \param[out] rho Array of nP density values, array must be valid before calling. void BlackoilPropertiesFromDeck::density(const int n, const double* A, + const int* cells, double* rho) const { const int np = numPhases(); - const double* sdens = pvt_.surfaceDensities(); // #pragma omp parallel for for (int i = 0; i < n; ++i) { + int cellIdx = cells?cells[i]:i; + int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx); + const double* sdens = pvt_.surfaceDensities(pvtRegionIdx); for (int phase = 0; phase < np; ++phase) { rho[np*i + phase] = 0.0; for (int comp = 0; comp < np; ++comp) { @@ -227,9 +242,10 @@ namespace Opm /// Densities of stock components at surface conditions. /// \return Array of P density values. - const double* BlackoilPropertiesFromDeck::surfaceDensity() const + const double* BlackoilPropertiesFromDeck::surfaceDensity(int cellIdx) const { - return pvt_.surfaceDensities(); + int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx); + return pvt_.surfaceDensities(pvtRegionIdx); } /// \param[in] n Number of data points. diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index e29194249..448c5cb56 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -97,6 +97,11 @@ namespace Opm /// \return N, the number of cells. virtual int numCells() const; + /// Return an array containing the PVT table index for each + /// grid cell + virtual const int* cellPvtRegionIndex() const + { return &cellPvtRegionIdx_[0]; } + /// \return Array of N porosity values. virtual const double* porosity() const; @@ -152,14 +157,16 @@ namespace Opm /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices /// are assumed to be in Fortran order, and are typically the result /// of a call to the method matrix(). + /// \param[in] cells The index of the grid cell of each data point. /// \param[out] rho Array of nP density values, array must be valid before calling. virtual void density(const int n, const double* A, + const int* cells, double* rho) const; /// Densities of stock components at surface conditions. /// \return Array of P density values. - virtual const double* surfaceDensity() const; + virtual const double* surfaceDensity(int cellIdx = 0) const; /// \param[in] n Number of data points. /// \param[in] s Array of nP saturation values. @@ -206,6 +213,13 @@ namespace Opm double* smax) const; private: + int getTableIndex_(const int* pvtTableIdx, int cellIdx) const + { + if (!pvtTableIdx) + return 0; + return pvtTableIdx[cellIdx]; + } + template void init(Opm::DeckConstPtr deck, int number_of_cells, @@ -224,6 +238,7 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock); RockFromDeck rock_; + std::vector cellPvtRegionIdx_; BlackoilPvtProperties pvt_; std::unique_ptr satprops_; mutable std::vector B_; diff --git a/opm/core/props/BlackoilPropertiesInterface.hpp b/opm/core/props/BlackoilPropertiesInterface.hpp index 80e3ef7bb..08edc8ae4 100644 --- a/opm/core/props/BlackoilPropertiesInterface.hpp +++ b/opm/core/props/BlackoilPropertiesInterface.hpp @@ -47,6 +47,10 @@ namespace Opm /// \return N, the number of cells. virtual int numCells() const = 0; + /// Return an array containing the PVT table index for each + /// grid cell + virtual const int* cellPvtRegionIndex() const = 0; + /// \return Array of N porosity values. virtual const double* porosity() const = 0; @@ -102,14 +106,16 @@ namespace Opm /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices /// are assumed to be in Fortran order, and are typically the result /// of a call to the method matrix(). + /// \param[in] cells The index of the grid cell of each data point. /// \param[out] rho Array of nP density values, array must be valid before calling. virtual void density(const int n, const double* A, + const int* cells, double* rho) const = 0; /// Densities of stock components at surface conditions. /// \return Array of P density values. - virtual const double* surfaceDensity() const = 0; + virtual const double* surfaceDensity(int regionIdx = 0) const = 0; /// \param[in] n Number of data points. /// \param[in] s Array of nP saturation values. diff --git a/opm/core/props/pvt/PvtPropertiesBasic.hpp b/opm/core/props/pvt/PvtPropertiesBasic.hpp index e9198ff97..46add6221 100644 --- a/opm/core/props/pvt/PvtPropertiesBasic.hpp +++ b/opm/core/props/pvt/PvtPropertiesBasic.hpp @@ -95,6 +95,8 @@ namespace Opm double* output_dRdp) const; private: + // The PVT properties. We need to store one value per PVT + // region. std::vector density_; std::vector viscosity_; std::vector formation_volume_factor_; diff --git a/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp b/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp index e23839e91..271763518 100644 --- a/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp +++ b/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp @@ -33,9 +33,9 @@ namespace Opm { } - void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck ) + void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck) { - // If we need multiple regions, this class and the Pvt* classes must change. + // So far, this class only supports a single PVT region. TODO? int region_number = 0; PhaseUsage phase_usage = phaseUsageFromDeck(deck); diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index dd90e3689..e97d5077f 100644 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ b/opm/core/simulator/EquilibrationHelpers.hpp @@ -139,7 +139,7 @@ namespace Opm props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp); std::vector rho(np, 0.0); - props_.density(1, &A[0], &rho[0]); + props_.density(1, &A[0], &c_[0], &rho[0]); return rho; } diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index e5f8c8ec2..a4ec1ee4f 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -369,8 +369,6 @@ namespace Opm const UnstructuredGrid& G , const double grav) { - typedef Miscibility::NoMixing NoMix; - for (typename RMap::RegionId r = 0, nr = reg.numRegions(); r < nr; ++r) diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index 6b1885ea1..70a6248dd 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -184,7 +184,7 @@ namespace Opm double A[4] = { 0.0 }; props_.matrix(1, &pressure, surfvol[phase], cells, A, 0); double rho[2] = { 0.0 }; - props_.density(1, A, rho); + props_.density(1, A, cells, rho); return rho[phase]; } }; diff --git a/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp index 4370a8ef5..82f8ea800 100644 --- a/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp +++ b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp @@ -418,7 +418,7 @@ namespace Opm assert(np == 2); const int dim = grid_.dimensions; density_.resize(nc*np); - props_.density(grid_.number_of_cells, &A_[0], &density_[0]); + props_.density(grid_.number_of_cells, &A_[0], /*cellIndices=*/NULL, &density_[0]); std::fill(gravflux_.begin(), gravflux_.end(), 0.0); for (int f = 0; f < nf; ++f) { const int* c = &grid_.face_cells[2*f]; diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 2454eb73a..8a0ebac8f 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -285,7 +285,8 @@ namespace Opm /// @brief Computes saturation from surface volume void computeSaturation(const BlackoilPropertiesInterface& props, - BlackoilState& state){ + BlackoilState& state) + { const int np = props.numPhases(); const int nc = props.numCells();