diff --git a/opm/autodiff/SolventPropsAdFromDeck.cpp b/opm/autodiff/SolventPropsAdFromDeck.cpp index 7bcb2f824..f4bd097e5 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.cpp +++ b/opm/autodiff/SolventPropsAdFromDeck.cpp @@ -44,6 +44,7 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, // retrieve the cell specific PVT table index from the deck // and using the grid... extractPvtTableIndex(cellPvtRegionIdx_, eclState, number_of_cells, global_cell); + extractTableIndex("SATNUM", cellSatNumRegionIdx_, eclState, number_of_cells, global_cell); // surface densities if (deck->hasKeyword("SDENSITY")) { @@ -84,7 +85,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, inverseBmu[i] = 1.0 / (b[i] * visc[i]); } - b_[regionIdx] = NonuniformTableLinear(press, inverseB); viscosity_[regionIdx] = NonuniformTableLinear(press, visc); inverseBmu_[regionIdx] = NonuniformTableLinear(press, inverseBmu); @@ -99,9 +99,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, int numRegions = ssfnTables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table saturation function supported for SSFN"); - } // resize the attributes of the object krg_.resize(numRegions); krs_.resize(numRegions); @@ -123,15 +120,18 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, if (deck->hasKeyword("MISCIBLE") ) { + + + // retrieve the cell specific Misc table index from the deck + // and using the grid... + extractTableIndex("MISCNUM", cellMiscRegionIdx_, eclState, number_of_cells, global_cell); + // misicible hydrocabon relative permeability wrt water const TableContainer& sof2Tables = tables->getSof2Tables(); if (!sof2Tables.empty()) { int numRegions = sof2Tables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table saturation function supported for SOF2"); - } // resize the attributes of the object krn_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { @@ -154,9 +154,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, int numRegions = miscTables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table miscibility function supported for MISC"); - } // resize the attributes of the object misc_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { @@ -180,9 +177,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, int numRegions = msfnTables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table saturation function supported for MSFN"); - } // resize the attributes of the object mkrsg_.resize(numRegions); mkro_.resize(numRegions); @@ -206,9 +200,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, int numRegions = sorwmisTables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table miscibility function supported for SORWMIS"); - } // resize the attributes of the object sorwmis_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { @@ -227,9 +218,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, int numRegions = sgcwmisTables.size(); - if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table miscibility function supported for SGCWMIS"); - } // resize the attributes of the object sgcwmis_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { @@ -244,12 +232,13 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, } if (deck->hasKeyword("TLMIXPAR")) { - const auto tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(0); - std::vector mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData(); - const int numRegions = mix_params_viscosity.size(); + const int numRegions = deck->getKeyword("TLMIXPAR")->size(); if (numRegions > 1) { OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR."); } + const auto tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(0); + std::vector mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData(); + mix_param_viscosity_ = mix_params_viscosity[0]; std::vector mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData(); @@ -282,7 +271,7 @@ ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg, V dmudp(n); for (int i = 0; i < n; ++i) { const double& pg_i = pg.value()[i]; - int regionIdx = cellPvtRegionIdx_[i]; + int regionIdx = cellPvtRegionIdx_[cells[i]]; double tempInvB = b_[regionIdx](pg_i); double tempInvBmu = inverseBmu_[regionIdx](pg_i); mu[i] = tempInvB / tempInvBmu; @@ -302,35 +291,35 @@ ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg, ADB SolventPropsAdFromDeck::bSolvent(const ADB& pg, const Cells& cells) const { - return SolventPropsAdFromDeck::makeADBfromTables(pg, cells, b_); + return SolventPropsAdFromDeck::makeADBfromTables(pg, cells, cellPvtRegionIdx_, b_); } ADB SolventPropsAdFromDeck::gasRelPermMultiplier(const ADB& solventFraction, const Cells& cells) const { - return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krg_); + return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellSatNumRegionIdx_, krg_); } ADB SolventPropsAdFromDeck::solventRelPermMultiplier(const ADB& solventFraction, const Cells& cells) const { - return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krs_); + return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellSatNumRegionIdx_, krs_); } ADB SolventPropsAdFromDeck::misicibleHydrocarbonWaterRelPerm(const ADB& Sn, const Cells& cells) const { - return SolventPropsAdFromDeck::makeADBfromTables(Sn, cells, krn_); + return SolventPropsAdFromDeck::makeADBfromTables(Sn, cells, cellSatNumRegionIdx_, krn_); } ADB SolventPropsAdFromDeck::miscibleSolventGasRelPermMultiplier(const ADB& Ssg, const Cells& cells) const { if (mkrsg_.size() > 0) { - return SolventPropsAdFromDeck::makeADBfromTables(Ssg, cells, mkrsg_); + return SolventPropsAdFromDeck::makeADBfromTables(Ssg, cells, cellSatNumRegionIdx_, mkrsg_); } // trivial function if not specified return Ssg; @@ -340,7 +329,7 @@ ADB SolventPropsAdFromDeck::miscibleOilRelPermMultiplier(const ADB& So, const Cells& cells) const { if (mkro_.size() > 0) { - return SolventPropsAdFromDeck::makeADBfromTables(So, cells, mkro_); + return SolventPropsAdFromDeck::makeADBfromTables(So, cells, cellSatNumRegionIdx_, mkro_); } // trivial function if not specified return So; @@ -350,14 +339,14 @@ ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction, const Cells& cells) const { - return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, misc_); + return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellMiscRegionIdx_, misc_); } ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw, const Cells& cells) const { if (sgcwmis_.size()>0) { - return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sgcwmis_); + return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, cellMiscRegionIdx_, sgcwmis_); } // return zeros if not specified return ADB::constant(V::Zero(Sw.size())); @@ -367,7 +356,7 @@ ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw, const Cells& cells) const { if (sorwmis_.size()>0) { - return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sorwmis_); + return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, cellMiscRegionIdx_, sorwmis_); } // return zeros if not specified return ADB::constant(V::Zero(Sw.size())); @@ -375,6 +364,7 @@ ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD, const Cells& cells, + const std::vector& regionIdx, const std::vector>& tables) const { const int n = cells.size(); assert(X_AD.value().size() == n); @@ -382,9 +372,8 @@ ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD, V dx(n); for (int i = 0; i < n; ++i) { const double& X_i = X_AD.value()[i]; - int regionIdx = 0; // TODO add mapping from cells to sat function table - x[i] = tables[regionIdx](X_i); - dx[i] = tables[regionIdx].derivative(X_i); + x[i] = tables[regionIdx[cells[i]]](X_i); + dx[i] = tables[regionIdx[cells[i]]].derivative(X_i); } ADB::M dx_diag(dx.matrix().asDiagonal()); @@ -402,7 +391,7 @@ V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const { const int n = cells.size(); V density(n); for (int i = 0; i < n; ++i) { - int regionIdx = cellPvtRegionIdx_[i]; + int regionIdx = cellPvtRegionIdx_[cells[i]]; density[i] = solvent_surface_densities_[regionIdx]; } return density; @@ -416,5 +405,24 @@ double SolventPropsAdFromDeck::mixingParameterDensity() const { return mix_param_density_; } +void SolventPropsAdFromDeck::extractTableIndex(const std::string& keyword, + std::vector &tableIdx, + Opm::EclipseStateConstPtr eclState, + size_t numCompressed, + const int *compressedToCartesianCellIdx) const +{ + //Get the Region data + const std::vector& regionData = eclState->getIntGridProperty(keyword)->getData(); + // Convert this into an array of compressed cells + // Eclipse uses Fortran-style indices which start at 1 + // instead of 0, we subtract 1. + tableIdx.resize(numCompressed); + for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) { + size_t cartesianCellIdx = compressedToCartesianCellIdx ? compressedToCartesianCellIdx[cellIdx]:cellIdx; + assert(cartesianCellIdx < regionData.size()); + tableIdx[cellIdx] = regionData[cartesianCellIdx] - 1; + } +} + } //namespace OPM diff --git a/opm/autodiff/SolventPropsAdFromDeck.hpp b/opm/autodiff/SolventPropsAdFromDeck.hpp index 5cc4d0dab..c2e439900 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.hpp +++ b/opm/autodiff/SolventPropsAdFromDeck.hpp @@ -139,14 +139,30 @@ private: /// Makes ADB from table values /// \param[in] X Array of n table lookup values. /// \param[in] cells Array of n cell indices to be associated with the fraction values. - /// \param[in] tables Vector of tables, one for each PVT region. + /// \param[in] tables Vector of tables, one for each PVT region. /// \return Array of n solvent density values. ADB makeADBfromTables(const ADB& X, const Cells& cells, + const std::vector& regionIdx, const std::vector>& tables) const; + /// Helper function to create an array containing the + /// table index of for each compressed cell from an Eclipse deck. + /// \param[in] keyword eclKeyword specifying region (SATNUM etc. ) + /// \param[in/out] tableIdx table index for each compressed cell + /// \param[in] eclState eclState from opm-parser + /// \param[in] numCompressed number of compressed cells + /// \param[in] compressedToCartesianCellIdx cartesianCellIdx for each cell in the grid + void extractTableIndex(const std::string& keyword, + std::vector &tableIdx, + Opm::EclipseStateConstPtr eclState, + size_t numCompressed, + const int *compressedToCartesianCellIdx) const; + // The PVT region which is to be used for each cell std::vector cellPvtRegionIdx_; + std::vector cellMiscRegionIdx_; + std::vector cellSatNumRegionIdx_; std::vector > b_; std::vector > viscosity_; std::vector > inverseBmu_;