OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #583 from totto82/solventSupportRegions

Browse Source 
Support regions in the solvent model
This commit is contained in:
Atgeirr Flø Rasmussen 2016-02-18 15:43:27 +01:00
commit 569df33ab2
4 changed files with 150 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
opm/autodiff/AutoDiffBlock.hpp
View File

@ -623,6 +623,36 @@ namespace Opm
return rhs * lhs; // Commutative operation.
}
/**
* @brief Computes the value of base raised to the power of exp elementwise
*
* @param base The AD forward block
* @param exp array of exponents
* @return The value of base raised to the power of exp elementwise
*/
template <typename Scalar>
AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
const typename AutoDiffBlock<Scalar>::V& exp)
{
const int num_elem = base.value().size();
typename AutoDiffBlock<Scalar>::V val (num_elem);
typename AutoDiffBlock<Scalar>::V derivative = exp;
assert(exp.size() == num_elem);
for (int i = 0; i < num_elem; ++i) {
val[i] = std::pow(base.value()[i], exp[i]);
derivative[i] *= std::pow(base.value()[i], exp[i] - 1.0);
}
const typename AutoDiffBlock<Scalar>::M derivative_diag(derivative.matrix().asDiagonal());
std::vector< typename AutoDiffBlock<Scalar>::M > jac (base.numBlocks());
for (int block = 0; block < base.numBlocks(); block++) {
fastSparseProduct(derivative_diag, base.derivative()[block], jac[block]);
}
return AutoDiffBlock<Scalar>::function( std::move(val), std::move(jac) );
}
/**
* @brief Computes the value of base raised to the power of exp
*
@ -635,7 +665,7 @@ namespace Opm
const double exp)
{
const typename AutoDiffBlock<Scalar>::V val = base.value().pow(exp);
const typename AutoDiffBlock<Scalar>::V derivative = exp * base.value().pow(exp-1.0);
const typename AutoDiffBlock<Scalar>::V derivative = exp * base.value().pow(exp - 1.0);
const typename AutoDiffBlock<Scalar>::M derivative_diag(derivative.matrix().asDiagonal());
std::vector< typename AutoDiffBlock<Scalar>::M > jac (base.numBlocks());

28
opm/autodiff/BlackoilSolventModel_impl.hpp
View File

@ -816,21 +816,21 @@ namespace Opm {
Selector<double> zero_selectorSsg(ssg_eff.value(), Selector<double>::Zero);
Selector<double> zero_selectorSn(sn_eff.value(), Selector<double>::Zero);
const ADB mu_s_pow = pow(mu_s,0.25);
const ADB mu_o_pow = pow(mu_o,0.25);
const ADB mu_g_pow = pow(mu_g,0.25);
const ADB mu_s_pow = pow(mu_s, 0.25);
const ADB mu_o_pow = pow(mu_o, 0.25);
const ADB mu_g_pow = pow(mu_g, 0.25);
const ADB mu_mos = zero_selectorSos.select(mu_o + mu_s, mu_o * mu_s / pow( ( (so_eff / sos_eff) * mu_s_pow) + ( (ss_eff / sos_eff) * mu_o_pow) , 4.0));
const ADB mu_msg = zero_selectorSsg.select(mu_g + mu_s , mu_g * mu_s / pow( ( (sg_eff / ssg_eff) * mu_s_pow) + ( (ss_eff / ssg_eff) * mu_g_pow) , 4.0));
const ADB mu_m = zero_selectorSn.select(mu_s + mu_o + mu_g, mu_o * mu_s * mu_g / pow( ( (so_eff / sn_eff) * mu_s_pow * mu_g_pow)
+ ( (ss_eff / sn_eff) * mu_o_pow * mu_g_pow) + ( (sg_eff / sn_eff) * mu_s_pow * mu_o_pow), 4.0));
// Mixing parameter for viscosity
const double mix_param_mu = solvent_props_.mixingParameterViscosity();
const V mix_param_mu = solvent_props_.mixingParameterViscosity(cells_);
// Update viscosities
viscosity[pu.phase_pos[ Oil ]] = pow(mu_o,1.0 - mix_param_mu) * pow(mu_mos,mix_param_mu);
viscosity[pu.phase_pos[ Gas ]] = pow(mu_g,1.0 - mix_param_mu) * pow(mu_msg,mix_param_mu);
viscosity[solvent_pos_] = pow(mu_s,1.0 - mix_param_mu) * pow(mu_m,mix_param_mu);
viscosity[pu.phase_pos[ Oil ]] = pow(mu_o,ones - mix_param_mu) * pow(mu_mos, mix_param_mu);
viscosity[pu.phase_pos[ Gas ]] = pow(mu_g,ones - mix_param_mu) * pow(mu_msg, mix_param_mu);
viscosity[solvent_pos_] = pow(mu_s,ones - mix_param_mu) * pow(mu_m, mix_param_mu);
// Density
ADB& rho_o = density[pu.phase_pos[ Oil ]];
@ -838,13 +838,13 @@ namespace Opm {
ADB& rho_s = density[solvent_pos_];
// mixing parameter for density
const double mix_param_rho = solvent_props_.mixingParameterDensity();
const V mix_param_rho = solvent_props_.mixingParameterDensity(cells_);
// compute effective viscosities for density calculations. These have to
// be recomputed as a different mixing parameter may be used.
const ADB mu_o_eff = pow(mu_o,1.0 - mix_param_rho) * pow(mu_mos,mix_param_rho);
const ADB mu_g_eff = pow(mu_g,1.0 - mix_param_rho) * pow(mu_msg,mix_param_rho);
const ADB mu_s_eff = pow(mu_s,1.0 - mix_param_rho) * pow(mu_m,mix_param_rho);
const ADB mu_o_eff = pow(mu_o,ones - mix_param_rho) * pow(mu_mos, mix_param_rho);
const ADB mu_g_eff = pow(mu_g,ones - mix_param_rho) * pow(mu_msg, mix_param_rho);
const ADB mu_s_eff = pow(mu_s,ones - mix_param_rho) * pow(mu_m, mix_param_rho);
const ADB sog_eff = so_eff + sg_eff;
const ADB sof = so_eff / sog_eff;
@ -852,9 +852,9 @@ namespace Opm {
// Effective densities
const ADB mu_sog_pow = mu_s_pow * ( (sgf * mu_o_pow) + (sof * mu_g_pow) );
const ADB mu_o_eff_pow = pow(mu_o_eff,0.25);
const ADB mu_g_eff_pow = pow(mu_g_eff,0.25);
const ADB mu_s_eff_pow = pow(mu_s_eff,0.25);
const ADB mu_o_eff_pow = pow(mu_o_eff, 0.25);
const ADB mu_g_eff_pow = pow(mu_g_eff, 0.25);
const ADB mu_s_eff_pow = pow(mu_s_eff, 0.25);
const ADB sfraction_oe = (mu_o_pow * (mu_o_eff_pow - mu_s_pow)) / (mu_o_eff_pow * (mu_o_pow - mu_s_pow));
const ADB sfraction_ge = (mu_g_pow * (mu_s_pow - mu_g_eff_pow)) / (mu_g_eff_pow * (mu_s_pow - mu_g_pow));
const ADB sfraction_se = (mu_sog_pow - ( mu_o_pow * mu_g_pow * mu_s_pow / mu_s_eff_pow) ) / ( mu_sog_pow - (mu_o_pow * mu_g_pow));

137
opm/autodiff/SolventPropsAdFromDeck.cpp
View File

@ -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", eclState, number_of_cells, global_cell, cellSatNumRegionIdx_);
// 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<double>(press, inverseB);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
inverseBmu_[regionIdx] = NonuniformTableLinear<double>(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", eclState, number_of_cells, global_cell, cellMiscRegionIdx_);
// 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,30 +232,27 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
}
if (deck->hasKeyword("TLMIXPAR")) {
const auto tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(0);
std::vector<double> mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData();
const int numRegions = mix_params_viscosity.size();
if (numRegions > 1) {
OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR.");
}
mix_param_viscosity_ = mix_params_viscosity[0];
const int numRegions = deck->getKeyword("TLMIXPAR")->size();
std::vector<double> mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData();
const int numDensityItems = mix_params_density.size();
if (numDensityItems == 0) {
mix_param_density_ = mix_param_viscosity_;
} else if (numDensityItems == 1) {
mix_param_density_ = mix_params_density[0];
} else {
OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR.");
// resize the attributes of the object
mix_param_viscosity_.resize(numRegions);
mix_param_density_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const auto& tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(regionIdx);
const auto& mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData();
mix_param_viscosity_[regionIdx] = mix_params_viscosity[0];
const auto& mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData();
const int numDensityItems = mix_params_density.size();
if (numDensityItems == 0) {
mix_param_density_[regionIdx] = mix_param_viscosity_[regionIdx];
} else if (numDensityItems == 1) {
mix_param_density_[regionIdx] = mix_params_density[0];
} else {
OPM_THROW(std::runtime_error, "Only one value can be entered for the TL parameter pr MISC region.");
}
}
} else {
mix_param_viscosity_ = 0.0;
mix_param_density_ = 0.0;
}
}
}
@ -282,7 +267,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 +287,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 +325,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 +335,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 +352,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 +360,7 @@ ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw
ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD,
const Cells& cells,
const std::vector<int>& regionIdx,
const std::vector<NonuniformTableLinear<double>>& tables) const {
const int n = cells.size();
assert(X_AD.value().size() == n);
@ -382,9 +368,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,18 +387,56 @@ 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;
}
double SolventPropsAdFromDeck::mixingParameterViscosity() const {
return mix_param_viscosity_;
V SolventPropsAdFromDeck::mixingParameterViscosity(const Cells& cells) const {
const int n = cells.size();
if (mix_param_viscosity_.size() > 0) {
V mix_param(n);
for (int i = 0; i < n; ++i) {
int regionIdx = cellMiscRegionIdx_[cells[i]];
mix_param[i] = mix_param_viscosity_[regionIdx];
}
return mix_param;
}
// return zeros if not specified
return V::Zero(n);
}
double SolventPropsAdFromDeck::mixingParameterDensity() const {
return mix_param_density_;
V SolventPropsAdFromDeck::mixingParameterDensity(const Cells& cells) const {
const int n = cells.size();
if (mix_param_viscosity_.size() > 0) {
V mix_param(n);
for (int i = 0; i < n; ++i) {
int regionIdx = cellMiscRegionIdx_[cells[i]];
mix_param[i] = mix_param_density_[regionIdx];
}
return mix_param;
}
// return zeros if not specified
return V::Zero(n);
}
void SolventPropsAdFromDeck::extractTableIndex(const std::string& keyword,
Opm::EclipseStateConstPtr eclState,
size_t numCompressed,
const int* compressedToCartesianCellIdx,
std::vector<int>& tableIdx) 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;
}
}

30
opm/autodiff/SolventPropsAdFromDeck.hpp
View File

@ -128,10 +128,14 @@ public:
V solventSurfaceDensity(const Cells& cells) const;
/// Todd-Longstaff mixing parameter for viscosity calculation
double mixingParameterViscosity() const;
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// return Array of n mixing paramters for viscosity calculation
V mixingParameterViscosity(const Cells& cells) const;
/// Todd-Longstaff mixing parameter for density calculation
double mixingParameterDensity() const;
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// return Array of n mixing paramters for density calculation
V mixingParameterDensity(const Cells& cells) const;
private:
@ -139,14 +143,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<int>& regionIdx,
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] eclState eclState from opm-parser
/// \param[in] numCompressed number of compressed cells
/// \param[in] compressedToCartesianCellIdx cartesianCellIdx for each cell in the grid
/// \param[out] tableIdx table index for each compressed cell
void extractTableIndex(const std::string& keyword,
Opm::EclipseStateConstPtr eclState,
size_t numCompressed,
const int* compressedToCartesianCellIdx,
std::vector<int>& tableIdx) const;
// The PVT region which is to be used for each cell
std::vector<int> cellPvtRegionIdx_;
std::vector<int> cellMiscRegionIdx_;
std::vector<int> cellSatNumRegionIdx_;
std::vector<NonuniformTableLinear<double> > b_;
std::vector<NonuniformTableLinear<double> > viscosity_;
std::vector<NonuniformTableLinear<double> > inverseBmu_;
@ -159,8 +179,8 @@ private:
std::vector<NonuniformTableLinear<double> > misc_;
std::vector<NonuniformTableLinear<double> > sorwmis_;
std::vector<NonuniformTableLinear<double> > sgcwmis_;
double mix_param_viscosity_;
double mix_param_density_;
std::vector<double> mix_param_viscosity_;
std::vector<double> mix_param_density_;
};
} // namespace OPM