Scaling of relperm functions - oil/water systems.

opm/core/eclipse/EclipseGridParser.cpp
opm/core/eclipse/EclipseGridParser.hpp
- New keywords: ENDSCALE SCALECRS SWCR SWL SWU SOWCR KRW KRWR KRO KRORW

opm/core/eclipse/SpecialEclipseFields.hpp
- Parsers for ENDSCALE and SCALECRS.

opm/core/fluid/BlackoilPropertiesFromDeck.cpp
- Consistency check: ENDSCALE implemented for SatFuncSimple only.

opm/core/fluid/SatFuncGwseg.hpp
opm/core/fluid/SatFuncSimple.hpp
opm/core/fluid/SatFuncStone2.hpp
- Accomodate "default" values for scalable parameters.
- For SatFuncGwseg and SatFuncStone2 the associated functionality not
  yet supported and the variables are dummies to satisfy the compiler.

opm/core/fluid/SatFuncSimple.cpp
- Initialisation for scalable parameters.
- Evaluation of relperms: Use (1-so) for evaluation of oil-relperms.
  (For scaled arguments sw and so do not necessarily add to one.)
- TODO: SatFuncGwseg.cpp and SatFuncStone2.cpp for oil-water systems.

opm/core/fluid/SaturationPropsFromDeck.hpp
- Struct to accomodate cell-wise scaling factors.
- Two flags indicating scaling and method.
- Methods for parameter initialisation and scaled relperm computation.

opm/core/fluid/SaturationPropsFromDeck_impl.hpp
- Initialize scaling options and relevant cell-wise scaling factors.
- Relperm evaluation modified for possible end point scaling.
This commit is contained in:
osae 2012-11-29 16:36:13 +01:00
parent 534f8f8bb6
commit 28a3a0770f
3 changed files with 303 additions and 2 deletions

View File

@ -50,6 +50,9 @@ namespace Opm
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 200);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
if (deck.hasField("ENDSCALE") && threephase_model != "simple") {
THROW("Sorry, end point scaling currently available for the 'simple' model only.");
}
if (sat_samples > 1) {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr

View File

@ -107,9 +107,28 @@ namespace Opm
std::vector<SatFuncSet> satfuncset_;
std::vector<int> cell_to_func_; // = SATNUM - 1
struct { // End point scaling parameters
std::vector<double> swl_;
std::vector<double> swcr_;
std::vector<double> swu_;
std::vector<double> sowcr_;
std::vector<double> krw_;
std::vector<double> krwr_;
std::vector<double> kro_;
std::vector<double> krorw_;
} eps_;
bool do_eps_; // ENDSCALE is active
bool do_3pt_; // SCALECRS: YES~true NO~false
typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
void initEPS(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const std::string& keyword,
std::vector<double>& scaleparam);
void relpermEPS(const double *s, const int cell, double *kr, double *dkrds= 0) const;
};

View File

@ -96,6 +96,61 @@ namespace Opm
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_, samples);
}
// Saturation table scaling
do_eps_ = false;
do_3pt_ = false;
if (deck.hasField("ENDSCALE")) {
if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) {
THROW("Currently endpoint-scaling limited to oil-water systems without gas.");
}
if (deck.getENDSCALE().dir_switch_ != std::string("NODIR")) {
THROW("SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'NODIR' accepted.");
}
if (deck.getENDSCALE().revers_switch_ != std::string("REVERS")) {
THROW("SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'REVERS' accepted.");
}
if (deck.hasField("SCALECRS")) {
if (deck.getSCALECRS().scalecrs_ == std::string("YES")) {
do_3pt_ = true;
}
}
do_eps_ = true;
initEPS(deck, grid, std::string("SWCR"), eps_.swcr_);
initEPS(deck, grid, std::string("SWL"), eps_.swl_);
initEPS(deck, grid, std::string("SWU"), eps_.swu_);
initEPS(deck, grid, std::string("SOWCR"), eps_.sowcr_);
initEPS(deck, grid, std::string("KRW"), eps_.krw_);
initEPS(deck, grid, std::string("KRWR"), eps_.krwr_);
initEPS(deck, grid, std::string("KRO"), eps_.kro_);
initEPS(deck, grid, std::string("KRORW"), eps_.krorw_);
/*
double ss[PhaseUsage::MaxNumPhases], kr[PhaseUsage::MaxNumPhases];
int oldP = std::cout.precision();
std::cout.precision(4);
for (unsigned int i=0; i<=100; ++i) {
ss[phase_usage_.phase_pos[Aqua]] = i*0.01;
ss[phase_usage_.phase_pos[Liquid]] = 1.0 - ss[phase_usage_.phase_pos[Aqua]];
endScaling(ss, 15, kr);
std::cout << std::showpoint
<< std::setw(10) << ss[phase_usage_.phase_pos[Aqua]]
<< std::setw(12) << kr[phase_usage_.phase_pos[Aqua]]
<< std::setw(10) << kr[phase_usage_.phase_pos[Liquid]];
endScaling(ss, 45, kr);
std::cout << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]]
<< std::setw(10) << kr[phase_usage_.phase_pos[Liquid]];
endScaling(ss, 75, kr);
std::cout << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]]
<< std::setw(10) << kr[phase_usage_.phase_pos[Liquid]];
endScaling(ss, 105, kr);
std::cout << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]]
<< std::setw(10) << kr[phase_usage_.phase_pos[Liquid]]
<< std::noshowpoint << std::endl;
}
std::cout.precision(oldP);
*/
}
}
@ -134,12 +189,20 @@ namespace Opm
if (dkrds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
if (do_eps_) {
relpermEPS(s + np*i, cells[i], kr + np*i, dkrds + np*np*i);
} else {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
}
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i);
if (do_eps_) {
relpermEPS(s + np*i, cells[i], kr + np*i);
} else {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i);
}
}
}
}
@ -214,7 +277,223 @@ namespace Opm
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
// Initialize saturation scaling parameter
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const std::string& keyword,
std::vector<double>& scaleparam)
{
if (deck.hasField(keyword)) {
// Active keyword assigned default values for each cell (in case of possible box-wise assignment)
scaleparam.resize(grid.number_of_cells);
int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua];
if (keyword == std::string("SWCR")) {
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).swcr_;
} else if (keyword == std::string("SWL")) {
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smin_[phase_pos_aqua];
} else if (keyword == std::string("SWU")) {
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smax_[phase_pos_aqua];
} else if (keyword == std::string("SOWCR")) {
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sowcr_;
} else if (keyword == std::string("KRW")) {
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwmax_;
} else if (keyword == std::string("KRWR")) {
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwr_;
} else if (keyword == std::string("KRO")) {
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).kromax_;
} else if (keyword == std::string("KRORW")) {
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorw_;
} else {
THROW("SaturationPropsFromDeck::initEndscale() -- unknown keyword: '" << keyword << "'");
}
// Keyword values from deck
std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
const int* gc = grid.global_cell;
const std::vector<double>& val = deck.getFloatingPointValue(keyword);
for (int c = 0; c < int(scaleparam.size()); ++c) {
const int deck_pos = (gc == NULL) ? c : gc[c];
scaleparam[c] = val[deck_pos];
}
}
}
// Saturation scaling
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::relpermEPS(const double *s, const int cell, double *kr, double *dkrds) const
{
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid];
double ss[PhaseUsage::MaxNumPhases];
if (do_3pt_) { // Three-point scaling
// Transforms for water saturation
if (eps_.swcr_.empty() && eps_.swu_.empty()) {
ss[wpos] = s[wpos];
} else {
double s_r = 1.0-funcForCell(cell).sowcr_;
double sr = eps_.sowcr_.empty() ? s_r : 1.0-eps_.sowcr_[cell];
if (s[wpos] <= sr) {
double sw_cr = funcForCell(cell).swcr_;
double swcr = eps_.swcr_.empty() ? sw_cr : eps_.swcr_[cell];
ss[wpos] = (s[wpos] <= swcr) ? sw_cr : sw_cr+(s[wpos]-swcr)*(s_r-sw_cr)/(sr-swcr);
} else {
double sw_max = funcForCell(cell).smax_[wpos];
double swmax = eps_.swu_.empty() ? sw_max : eps_.swu_[cell];
ss[wpos] = (s[wpos] >= swmax) ? sw_max : s_r+(s[wpos]-sr)*(sw_max-s_r)/(swmax-sr);
}
}
// Transforms for oil saturation
if (eps_.sowcr_.empty() && eps_.swl_.empty()) {
ss[opos] = s[opos];
} else {
double s_r = 1.0-funcForCell(cell).swcr_;
double sr = eps_.swcr_.empty() ? s_r : 1.0-eps_.swcr_[cell];
if (s[opos] <= sr) {
double sow_cr = funcForCell(cell).sowcr_;
double sowcr = eps_.sowcr_.empty() ? sow_cr : eps_.sowcr_[cell];
ss[opos] = (s[opos] <= sowcr) ? sow_cr : sow_cr+(s[opos]-sowcr)*(s_r-sow_cr)/(sr-sowcr);
} else {
double sow_max = funcForCell(cell).smax_[opos];
double sowmax = eps_.swl_.empty() ? sow_max : (1.0-eps_.swl_[cell]);
ss[opos] = (s[opos] >= sowmax) ? sow_max : s_r+(s[opos]-sr)*(sow_max-s_r)/(sowmax-sr);
}
}
} else { // Two-point scaling
// Transforms for water saturation
if (eps_.swcr_.empty() && eps_.swu_.empty()) {
ss[wpos] = s[wpos];
} else {
double sw_cr = funcForCell(cell).swcr_;
double swcr = eps_.swcr_.empty() ? sw_cr : eps_.swcr_[cell];
if (s[wpos] <= swcr) {
ss[wpos] = sw_cr;
} else {
double sw_max = funcForCell(cell).smax_[wpos];
double swmax = eps_.swu_.empty() ? sw_max : eps_.swu_[cell];
ss[wpos] = (s[wpos] >= swmax) ? sw_max : sw_cr + (s[wpos]-swcr)*(sw_max-sw_cr)/(swmax-swcr);
}
}
// Transforms for oil saturation
if (eps_.sowcr_.empty() && eps_.swl_.empty()) {
ss[opos] = s[opos];
} else {
double sow_cr = funcForCell(cell).sowcr_;
double socr = eps_.sowcr_.empty() ? sow_cr : eps_.sowcr_[cell];
if (s[opos] <= socr) {
ss[opos] = sow_cr;
} else {
double sow_max = funcForCell(cell).smax_[opos];
double sowmax = eps_.swl_.empty() ? sow_max : (1.0-eps_.swl_[cell]);
ss[opos] = (s[opos] >= sowmax) ? sow_max : sow_cr + (s[opos]-socr) *(sow_max-sow_cr)/(sowmax-socr);
}
}
}
// Evaluation of relperms
if (dkrds) {
THROW("Relperm derivatives not yet available in combination with end point scaling ...");
funcForCell(cell).evalKrDeriv(ss, kr, dkrds);
} else {
// Assume: sw_cr -> krw=0 sw_max -> krw=<max water relperm>
// sow_cr -> kro=0 sow_max -> kro=<max oil relperm>
funcForCell(cell).evalKr(ss, kr);
}
// Scaling of relperms values
// - Water
if (eps_.krw_.empty() && eps_.krwr_.empty()) { // No value scaling
} else if (eps_.krwr_.empty()) { // Two-point
kr[wpos] *= (eps_.krw_[cell]/funcForCell(cell).krwmax_);
} else {
double swcr = eps_.swcr_.empty() ? funcForCell(cell).swcr_ : eps_.swcr_[cell];
double swmax = eps_.swu_.empty() ? funcForCell(cell).smax_[wpos] : eps_.swu_[cell];
double sr;
if (do_3pt_) {
sr = eps_.sowcr_.empty() ? 1.0-funcForCell(cell).sowcr_ : 1.0-eps_.sowcr_[cell];
} else {
double sw_cr = funcForCell(cell).swcr_;
double sw_max = funcForCell(cell).smax_[wpos];
double s_r = 1.0-funcForCell(cell).sowcr_;
sr = swcr + (s_r-sw_cr)*(swmax-swcr)/(sw_max-sw_cr);
}
if (s[wpos] <= swcr) {
kr[wpos] = 0.0;
} else if (sr > swmax-1.0e-6) {
if (do_3pt_) { //Ignore krw and do two-point?
kr[wpos] *= eps_.krwr_[cell]/funcForCell(cell).krwr_;
} else if (!eps_.kro_.empty()){ //Ignore krwr and do two-point
kr[wpos] *= eps_.krw_[cell]/funcForCell(cell).krwmax_;
}
} else if (s[wpos] <= sr) {
kr[wpos] *= eps_.krwr_[cell]/funcForCell(cell).krwr_;
} else if (s[wpos] <= swmax) {
double krw_max = funcForCell(cell).krwmax_;
double krw = eps_.krw_.empty() ? krw_max : eps_.krw_[cell];
double krw_r = funcForCell(cell).krwr_;
double krwr = eps_.krwr_.empty() ? krw_r : eps_.krwr_[cell];
if (std::fabs(krw_max- krw_r) > 1.0e-6) {
kr[wpos] = krwr + (kr[wpos]-krw_r)*(krw-krwr)/(krw_max-krw_r);
} else {
kr[wpos] = krwr + (krw-krwr)*(s[wpos]-sr)/(swmax-sr);
}
} else {
kr[wpos] = eps_.krw_.empty() ? funcForCell(cell).krwmax_ : eps_.krw_[cell];
}
}
// - Oil
if (eps_.kro_.empty() && eps_.krorw_.empty()) { // No value scaling
} else if (eps_.krorw_.empty()) { // Two-point scaling
kr[opos] *= (eps_.kro_[cell]/funcForCell(cell).kromax_);
} else {
double sowcr = eps_.sowcr_.empty() ? funcForCell(cell).sowcr_ : eps_.sowcr_[cell];
double sowmax = eps_.swl_.empty() ? funcForCell(cell).smax_[opos] : 1.0-eps_.swl_[cell];
double sr;
if (do_3pt_) {
sr = eps_.swcr_.empty() ? 1.0-funcForCell(cell).swcr_ : 1.0-eps_.swcr_[cell];
} else {
double sow_cr = funcForCell(cell).sowcr_;
double sow_max = funcForCell(cell).smax_[opos];
double s_r = 1.0-funcForCell(cell).swcr_;
sr = sowcr + (s_r-sow_cr)*(sowmax-sowcr)/(sow_max-sow_cr);
}
if (s[opos] <= sowcr) {
kr[opos] = 0.0;
} else if (sr > sowmax-1.0e-6) {
if (do_3pt_) { //Ignore kro and do two-point?
kr[opos] *= eps_.krorw_[cell]/funcForCell(cell).krorw_;
} else if (!eps_.kro_.empty()){ //Ignore krowr and do two-point
kr[opos] *= eps_.kro_[cell]/funcForCell(cell).kromax_;
}
} else if (s[opos] <= sr) {
kr[opos] *= eps_.krorw_[cell]/funcForCell(cell).krorw_;
} else if (s[opos] <= sowmax) {
double kro_max = funcForCell(cell).kromax_;
double kro = eps_.kro_.empty() ? kro_max : eps_.kro_[cell];
double kro_rw = funcForCell(cell).krorw_;
double krorw = eps_.krorw_[cell];
if (std::fabs(kro_max- kro_rw) > 1.0e-6) {
kr[opos] = krorw + (kr[opos]- kro_rw)*(kro-krorw)/(kro_max- kro_rw);
} else {
kr[opos] = krorw + (kro-krorw)*(s[opos]-sr)/(sowmax-sr);
}
} else {
kr[opos] = eps_.kro_.empty() ? funcForCell(cell).kromax_ : eps_.kro_[cell];
}
}
}
} // namespace Opm