Support pvt and saturation regions in SolventPropsAdFromDeck

This commit is contained in:
Tor Harald Sandve 2016-02-15 10:07:47 +01:00
parent eb228835c2
commit 9328bd5af5
2 changed files with 61 additions and 37 deletions

View File

@ -44,6 +44,7 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
// retrieve the cell specific PVT table index from the deck // retrieve the cell specific PVT table index from the deck
// and using the grid... // and using the grid...
extractPvtTableIndex(cellPvtRegionIdx_, eclState, number_of_cells, global_cell); extractPvtTableIndex(cellPvtRegionIdx_, eclState, number_of_cells, global_cell);
extractTableIndex("SATNUM", cellSatNumRegionIdx_, eclState, number_of_cells, global_cell);
// surface densities // surface densities
if (deck->hasKeyword("SDENSITY")) { if (deck->hasKeyword("SDENSITY")) {
@ -84,7 +85,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
inverseBmu[i] = 1.0 / (b[i] * visc[i]); inverseBmu[i] = 1.0 / (b[i] * visc[i]);
} }
b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB); b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc); viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu); inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu);
@ -99,9 +99,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
int numRegions = ssfnTables.size(); 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 // resize the attributes of the object
krg_.resize(numRegions); krg_.resize(numRegions);
krs_.resize(numRegions); krs_.resize(numRegions);
@ -123,15 +120,18 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
if (deck->hasKeyword("MISCIBLE") ) { 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 // misicible hydrocabon relative permeability wrt water
const TableContainer& sof2Tables = tables->getSof2Tables(); const TableContainer& sof2Tables = tables->getSof2Tables();
if (!sof2Tables.empty()) { if (!sof2Tables.empty()) {
int numRegions = sof2Tables.size(); 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 // resize the attributes of the object
krn_.resize(numRegions); krn_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
@ -154,9 +154,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
int numRegions = miscTables.size(); 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 // resize the attributes of the object
misc_.resize(numRegions); misc_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
@ -180,9 +177,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
int numRegions = msfnTables.size(); 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 // resize the attributes of the object
mkrsg_.resize(numRegions); mkrsg_.resize(numRegions);
mkro_.resize(numRegions); mkro_.resize(numRegions);
@ -206,9 +200,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
int numRegions = sorwmisTables.size(); 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 // resize the attributes of the object
sorwmis_.resize(numRegions); sorwmis_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
@ -227,9 +218,6 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
int numRegions = sgcwmisTables.size(); 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 // resize the attributes of the object
sgcwmis_.resize(numRegions); sgcwmis_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
@ -244,12 +232,13 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
} }
if (deck->hasKeyword("TLMIXPAR")) { if (deck->hasKeyword("TLMIXPAR")) {
const auto tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(0); const int numRegions = deck->getKeyword("TLMIXPAR")->size();
std::vector<double> mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData();
const int numRegions = mix_params_viscosity.size();
if (numRegions > 1) { if (numRegions > 1) {
OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR."); OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR.");
} }
const auto tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(0);
std::vector<double> mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData();
mix_param_viscosity_ = mix_params_viscosity[0]; mix_param_viscosity_ = mix_params_viscosity[0];
std::vector<double> mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData(); std::vector<double> mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData();
@ -282,7 +271,7 @@ ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg,
V dmudp(n); V dmudp(n);
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
const double& pg_i = pg.value()[i]; const double& pg_i = pg.value()[i];
int regionIdx = cellPvtRegionIdx_[i]; int regionIdx = cellPvtRegionIdx_[cells[i]];
double tempInvB = b_[regionIdx](pg_i); double tempInvB = b_[regionIdx](pg_i);
double tempInvBmu = inverseBmu_[regionIdx](pg_i); double tempInvBmu = inverseBmu_[regionIdx](pg_i);
mu[i] = tempInvB / tempInvBmu; mu[i] = tempInvB / tempInvBmu;
@ -302,35 +291,35 @@ ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg,
ADB SolventPropsAdFromDeck::bSolvent(const ADB& pg, ADB SolventPropsAdFromDeck::bSolvent(const ADB& pg,
const Cells& cells) const const Cells& cells) const
{ {
return SolventPropsAdFromDeck::makeADBfromTables(pg, cells, b_); return SolventPropsAdFromDeck::makeADBfromTables(pg, cells, cellPvtRegionIdx_, b_);
} }
ADB SolventPropsAdFromDeck::gasRelPermMultiplier(const ADB& solventFraction, ADB SolventPropsAdFromDeck::gasRelPermMultiplier(const ADB& solventFraction,
const Cells& cells) const const Cells& cells) const
{ {
return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krg_); return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellSatNumRegionIdx_, krg_);
} }
ADB SolventPropsAdFromDeck::solventRelPermMultiplier(const ADB& solventFraction, ADB SolventPropsAdFromDeck::solventRelPermMultiplier(const ADB& solventFraction,
const Cells& cells) const const Cells& cells) const
{ {
return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krs_); return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellSatNumRegionIdx_, krs_);
} }
ADB SolventPropsAdFromDeck::misicibleHydrocarbonWaterRelPerm(const ADB& Sn, ADB SolventPropsAdFromDeck::misicibleHydrocarbonWaterRelPerm(const ADB& Sn,
const Cells& cells) const const Cells& cells) const
{ {
return SolventPropsAdFromDeck::makeADBfromTables(Sn, cells, krn_); return SolventPropsAdFromDeck::makeADBfromTables(Sn, cells, cellSatNumRegionIdx_, krn_);
} }
ADB SolventPropsAdFromDeck::miscibleSolventGasRelPermMultiplier(const ADB& Ssg, ADB SolventPropsAdFromDeck::miscibleSolventGasRelPermMultiplier(const ADB& Ssg,
const Cells& cells) const const Cells& cells) const
{ {
if (mkrsg_.size() > 0) { if (mkrsg_.size() > 0) {
return SolventPropsAdFromDeck::makeADBfromTables(Ssg, cells, mkrsg_); return SolventPropsAdFromDeck::makeADBfromTables(Ssg, cells, cellSatNumRegionIdx_, mkrsg_);
} }
// trivial function if not specified // trivial function if not specified
return Ssg; return Ssg;
@ -340,7 +329,7 @@ ADB SolventPropsAdFromDeck::miscibleOilRelPermMultiplier(const ADB& So,
const Cells& cells) const const Cells& cells) const
{ {
if (mkro_.size() > 0) { if (mkro_.size() > 0) {
return SolventPropsAdFromDeck::makeADBfromTables(So, cells, mkro_); return SolventPropsAdFromDeck::makeADBfromTables(So, cells, cellSatNumRegionIdx_, mkro_);
} }
// trivial function if not specified // trivial function if not specified
return So; return So;
@ -350,14 +339,14 @@ ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction,
const Cells& cells) const const Cells& cells) const
{ {
return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, misc_); return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellMiscRegionIdx_, misc_);
} }
ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw, ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw,
const Cells& cells) const { const Cells& cells) const {
if (sgcwmis_.size()>0) { if (sgcwmis_.size()>0) {
return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sgcwmis_); return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, cellMiscRegionIdx_, sgcwmis_);
} }
// return zeros if not specified // return zeros if not specified
return ADB::constant(V::Zero(Sw.size())); return ADB::constant(V::Zero(Sw.size()));
@ -367,7 +356,7 @@ ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw
ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw, ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw,
const Cells& cells) const { const Cells& cells) const {
if (sorwmis_.size()>0) { if (sorwmis_.size()>0) {
return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sorwmis_); return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, cellMiscRegionIdx_, sorwmis_);
} }
// return zeros if not specified // return zeros if not specified
return ADB::constant(V::Zero(Sw.size())); return ADB::constant(V::Zero(Sw.size()));
@ -375,6 +364,7 @@ ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw
ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD, ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD,
const Cells& cells, const Cells& cells,
const std::vector<int>& regionIdx,
const std::vector<NonuniformTableLinear<double>>& tables) const { const std::vector<NonuniformTableLinear<double>>& tables) const {
const int n = cells.size(); const int n = cells.size();
assert(X_AD.value().size() == n); assert(X_AD.value().size() == n);
@ -382,9 +372,8 @@ ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD,
V dx(n); V dx(n);
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
const double& X_i = X_AD.value()[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[cells[i]]](X_i);
x[i] = tables[regionIdx](X_i); dx[i] = tables[regionIdx[cells[i]]].derivative(X_i);
dx[i] = tables[regionIdx].derivative(X_i);
} }
ADB::M dx_diag(dx.matrix().asDiagonal()); ADB::M dx_diag(dx.matrix().asDiagonal());
@ -402,7 +391,7 @@ V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const {
const int n = cells.size(); const int n = cells.size();
V density(n); V density(n);
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
int regionIdx = cellPvtRegionIdx_[i]; int regionIdx = cellPvtRegionIdx_[cells[i]];
density[i] = solvent_surface_densities_[regionIdx]; density[i] = solvent_surface_densities_[regionIdx];
} }
return density; return density;
@ -416,5 +405,24 @@ double SolventPropsAdFromDeck::mixingParameterDensity() const {
return mix_param_density_; return mix_param_density_;
} }
void SolventPropsAdFromDeck::extractTableIndex(const std::string& keyword,
std::vector<int> &tableIdx,
Opm::EclipseStateConstPtr eclState,
size_t numCompressed,
const int *compressedToCartesianCellIdx) const
{
//Get the Region data
const std::vector<int>& 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 } //namespace OPM

View File

@ -139,14 +139,30 @@ private:
/// Makes ADB from table values /// Makes ADB from table values
/// \param[in] X Array of n table lookup 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] 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. /// \return Array of n solvent density values.
ADB makeADBfromTables(const ADB& X, ADB makeADBfromTables(const ADB& X,
const Cells& cells, const Cells& cells,
const std::vector<int>& regionIdx,
const std::vector<NonuniformTableLinear<double>>& tables) const; const std::vector<NonuniformTableLinear<double>>& 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<int> &tableIdx,
Opm::EclipseStateConstPtr eclState,
size_t numCompressed,
const int *compressedToCartesianCellIdx) const;
// The PVT region which is to be used for each cell // The PVT region which is to be used for each cell
std::vector<int> cellPvtRegionIdx_; std::vector<int> cellPvtRegionIdx_;
std::vector<int> cellMiscRegionIdx_;
std::vector<int> cellSatNumRegionIdx_;
std::vector<NonuniformTableLinear<double> > b_; std::vector<NonuniformTableLinear<double> > b_;
std::vector<NonuniformTableLinear<double> > viscosity_; std::vector<NonuniformTableLinear<double> > viscosity_;
std::vector<NonuniformTableLinear<double> > inverseBmu_; std::vector<NonuniformTableLinear<double> > inverseBmu_;