OilfieldToolsNet/opm-core
4
0 Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

SaturationPropsFromDeck: make the jump to fluid states

Browse Source 
this means the following changes:

- the "SatFuncGwseg" class is converted
- for now, Gwseg is the only saturation function supported by
  SaturationPropsFromDeck. (will be changed in later commits.)
- the funcForCell() method of SaturationPropsFromDeck is removed as it
  just occludes things
This commit is contained in:
Andreas Lauser committed 2015-06-29 14:49:20 +02:00
1 parent 4b5de05d15
commit 63075249a9
3 files changed
+195 -161
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

opm/core/props/satfunc/SatFuncGwseg.hpp
+76 -59
View File
@@ -27,34 +27,42 @@ namespace Opm
class SatFuncGwseg : public SatFuncBase<TableType> class SatFuncGwseg : public SatFuncBase<TableType>
{ {
public: public:
void evalKr(const double* s, double* kr) const; template <class FluidState>
void evalKrDeriv(const double* s, double* kr, double* dkrds) const; void evalKr(const FluidState& fluidState, double* kr) const;
void evalPc(const double* s, double* pc) const; template <class FluidState>
void evalPcDeriv(const double* s, double* pc, double* dpcds) const; void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const;
template <class FluidState>
void evalKr(const double* s, double* kr, const EPSTransforms* epst) const; void evalPc(const FluidState& fluidState, double* pc) const;
void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const; template <class FluidState>
void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const; void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const;
void evalPc(const double* s, double* pc, const EPSTransforms* epst) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const;
private:
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst) const;
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const;
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const;
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const;
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc, const EPSTransforms* epst) const;
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds, const EPSTransforms* epst) const;
}; };
typedef SatFuncGwseg<UniformTableLinear<double> > SatFuncGwsegUniform; typedef SatFuncGwseg<UniformTableLinear<double> > SatFuncGwsegUniform;
typedef SatFuncGwseg<NonuniformTableLinear<double> > SatFuncGwsegNonuniform; typedef SatFuncGwseg<NonuniformTableLinear<double> > SatFuncGwsegNonuniform;
template<class TableType> template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, double* kr) const template <class FluidState>
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr) const
{ {
if (this->phase_usage.num_phases == 3) { if (this->phase_usage.num_phases == 3) {
// Relative permeability model based on segregation of water // Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones. // and gas, with oil present in both water and gas zones.
double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
const double sw = std::max(s[BlackoilPhases::Aqua], swco); const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = s[BlackoilPhases::Vapour]; const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-5; const double eps = 1e-5;
swco = std::min(swco,sw-eps); swco = std::min(swco,sw-eps);
@@ -77,7 +85,7 @@ namespace Opm
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos]; double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw); double krw = this->krw_(sw);
double krow = this->krow_(sw); double krow = this->krow_(sw);
kr[wpos] = krw; kr[wpos] = krw;
@@ -86,7 +94,7 @@ namespace Opm
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos]; double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg); double krg = this->krg_(sg);
double krog = this->krog_(sg); double krog = this->krog_(sg);
kr[gpos] = krg; kr[gpos] = krg;
@@ -96,7 +104,8 @@ namespace Opm
template<class TableType> template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, double* kr, const EPSTransforms* epst) const template <class FluidState>
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst) const
{ {
if (this->phase_usage.num_phases == 3) { if (this->phase_usage.num_phases == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
@@ -108,8 +117,8 @@ namespace Opm
// TODO Also consider connate gas ... // TODO Also consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin; double swco = epst->wat.smin;
const double sw = std::max(s[BlackoilPhases::Aqua], swco); const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = s[BlackoilPhases::Vapour]; const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6; const double eps = 1e-6;
swco = std::min(swco,sw-eps); swco = std::min(swco,sw-eps);
const double ssw = sg + sw; const double ssw = sg + sw;
@@ -141,7 +150,7 @@ namespace Opm
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos]; double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw); double krw = this->krw_(sw);
double krow = this->krow_(sw); double krow = this->krow_(sw);
kr[wpos] = krw; kr[wpos] = krw;
@@ -150,7 +159,7 @@ namespace Opm
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos]; double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg); double krg = this->krg_(sg);
double krog = this->krog_(sg); double krog = this->krog_(sg);
kr[gpos] = krg; kr[gpos] = krg;
@@ -159,7 +168,8 @@ namespace Opm
} }
template<class TableType> template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const template <class FluidState>
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{ {
if (this->phase_usage.num_phases == 3) { if (this->phase_usage.num_phases == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
@@ -171,8 +181,8 @@ namespace Opm
// TODO Consider connate gas ... // TODO Consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin; double swco = epst->wat.smin;
const double sw = std::max(s[BlackoilPhases::Aqua], swco); const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = s[BlackoilPhases::Vapour]; const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6; const double eps = 1e-6;
swco = std::min(swco,sw-eps); swco = std::min(swco,sw-eps);
const double ssw = sg + sw; const double ssw = sg + sw;
@@ -233,7 +243,7 @@ namespace Opm
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos]; double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw); double krw = this->krw_(sw);
double krow = this->krow_(sw); double krow = this->krow_(sw);
kr[wpos] = krw; kr[wpos] = krw;
@@ -242,7 +252,7 @@ namespace Opm
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos]; double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg); double krg = this->krg_(sg);
double krog = this->krog_(sg); double krog = this->krog_(sg);
kr[gpos] = krg; kr[gpos] = krg;
@@ -251,7 +261,8 @@ namespace Opm
} }
template<class TableType> template<class TableType>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds) const template<class FluidState>
void SatFuncGwseg<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
{ {
const int np = this->phase_usage.num_phases; const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0); std::fill(dkrds, dkrds + np*np, 0.0);
@@ -261,8 +272,8 @@ namespace Opm
// and gas, with oil present in both water and gas zones. // and gas, with oil present in both water and gas zones.
double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
const double sw = std::max(s[BlackoilPhases::Aqua], swco); const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = s[BlackoilPhases::Vapour]; const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-5; const double eps = 1e-5;
swco = std::min(swco,sw-eps); swco = std::min(swco,sw-eps);
@@ -297,7 +308,7 @@ namespace Opm
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos]; double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw); double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw); double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw); double krow = this->krow_(sw);
@@ -310,7 +321,7 @@ namespace Opm
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos]; double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg); double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg); double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg); double krog = this->krog_(sg);
@@ -324,7 +335,8 @@ namespace Opm
} }
template<class TableType> template<class TableType>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const template<class FluidState>
void SatFuncGwseg<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const
{ {
const int np = this->phase_usage.num_phases; const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0); std::fill(dkrds, dkrds + np*np, 0.0);
@@ -339,8 +351,8 @@ namespace Opm
// TODO Also consider connate gas ... // TODO Also consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin; double swco = epst->wat.smin;
const double sw = std::max(s[BlackoilPhases::Aqua], swco); const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = s[BlackoilPhases::Vapour]; const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6; const double eps = 1e-6;
swco = std::min(swco,sw-eps); swco = std::min(swco,sw-eps);
const double ssw = sg + sw; const double ssw = sg + sw;
@@ -386,7 +398,7 @@ namespace Opm
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos]; double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw); double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw); double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw); double krow = this->krow_(sw);
@@ -399,7 +411,7 @@ namespace Opm
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos]; double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg); double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg); double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg); double krog = this->krog_(sg);
@@ -412,7 +424,8 @@ namespace Opm
} }
template<class TableType> template<class TableType>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const template<class FluidState>
void SatFuncGwseg<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{ {
const int np = this->phase_usage.num_phases; const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0); std::fill(dkrds, dkrds + np*np, 0.0);
@@ -427,8 +440,8 @@ namespace Opm
// TODO Also consider connate gas ... // TODO Also consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin; double swco = epst->wat.smin;
const double sw = std::max(s[BlackoilPhases::Aqua], swco); const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = s[BlackoilPhases::Vapour]; const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6; const double eps = 1e-6;
swco = std::min(swco,sw-eps); swco = std::min(swco,sw-eps);
const double ssw = sg + sw; const double ssw = sg + sw;
@@ -506,7 +519,7 @@ namespace Opm
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos]; double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw); double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw); double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw); double krow = this->krow_(sw);
@@ -519,7 +532,7 @@ namespace Opm
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos]; double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg); double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg); double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg); double krog = this->krog_(sg);
@@ -532,38 +545,41 @@ namespace Opm
} }
template<class TableType> template<class TableType>
void SatFuncGwseg<TableType>::evalPc(const double* s, double* pc) const template <class FluidState>
void SatFuncGwseg<TableType>::evalPc(const FluidState& fluidState, double* pc) const
{ {
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]); pc[pos] = this->pcow_(fluidState.saturation(pos));
} }
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]); pc[pos] = this->pcog_(fluidState.saturation(pos));
} }
} }
template<class TableType> template<class TableType>
void SatFuncGwseg<TableType>::evalPc(const double* s, double* pc, const EPSTransforms* epst) const template <class FluidState>
void SatFuncGwseg<TableType>::evalPc(const FluidState& fluidState, double* pc, const EPSTransforms* epst) const
{ {
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
double _sw = epst->wat.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]); double _sw = epst->wat.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->wat.pcFactor*this->pcow_(_sw); pc[pos] = epst->wat.pcFactor*this->pcow_(_sw);
} }
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
double _sg = epst->gas.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]); double _sg = epst->gas.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->gas.pcFactor*this->pcog_(_sg); pc[pos] = epst->gas.pcFactor*this->pcog_(_sg);
} }
} }
template<class TableType> template<class TableType>
void SatFuncGwseg<TableType>::evalPcDeriv(const double* s, double* pc, double* dpcds) const template <class FluidState>
void SatFuncGwseg<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
{ {
// The problem of determining three-phase capillary pressures // The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase // is very hard experimentally, usually one extends two-phase
@@ -576,19 +592,20 @@ namespace Opm
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]); pc[pos] = this->pcow_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcow_.derivative(s[pos]); dpcds[np*pos + pos] = this->pcow_.derivative(fluidState.saturation(pos));
} }
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]); pc[pos] = this->pcog_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]); dpcds[np*pos + pos] = this->pcog_.derivative(fluidState.saturation(pos));
} }
} }
template<class TableType> template<class TableType>
void SatFuncGwseg<TableType>::evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const template <class FluidState>
void SatFuncGwseg<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds, const EPSTransforms* epst) const
{ {
// The problem of determining three-phase capillary pressures // The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase // is very hard experimentally, usually one extends two-phase
@@ -601,16 +618,16 @@ namespace Opm
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
double _sw = epst->wat.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]); double _sw = epst->wat.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->wat.pcFactor*this->pcow_(_sw); pc[pos] = epst->wat.pcFactor*this->pcow_(_sw);
double _dsdsw = epst->wat.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]); double _dsdsw = epst->wat.scaleSatDerivPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
dpcds[np*pos + pos] = epst->wat.pcFactor*_dsdsw*this->pcow_.derivative(_sw); dpcds[np*pos + pos] = epst->wat.pcFactor*_dsdsw*this->pcow_.derivative(_sw);
} }
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
double _sg = epst->gas.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]); double _sg = epst->gas.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->gas.pcFactor*this->pcog_(_sg); pc[pos] = epst->gas.pcFactor*this->pcog_(_sg);
double _dsdsg = epst->gas.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]); double _dsdsg = epst->gas.scaleSatDerivPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
dpcds[np*pos + pos] = epst->gas.pcFactor*_dsdsg*this->pcog_.derivative(_sg); dpcds[np*pos + pos] = epst->gas.pcFactor*_dsdsg*this->pcog_.derivative(_sg);
} }
} }
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
+9 -7
View File
@@ -23,8 +23,6 @@
#include <opm/core/props/satfunc/SaturationPropsInterface.hpp> #include <opm/core/props/satfunc/SaturationPropsInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/BlackoilPhases.hpp> #include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SatFuncStone2.hpp>
#include <opm/core/props/satfunc/SatFuncSimple.hpp>
#include <opm/core/props/satfunc/SatFuncGwseg.hpp> #include <opm/core/props/satfunc/SatFuncGwseg.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp> #include <opm/parser/eclipse/Deck/Deck.hpp>
@@ -129,10 +127,17 @@ namespace Opm
double & swat); double & swat);
private: private:
typedef SatFuncGwsegNonuniform SatFuncSet; // internal helper method for satRange()
template <class SaturationFunction>
void satRange_(const SaturationFunction& satFunc,
const int cellIdx,
const int* cells,
double* smin,
double* smax) const;
PhaseUsage phase_usage_; PhaseUsage phase_usage_;
std::vector<SatFuncSet> satfuncset_; typedef Opm::SatFuncGwseg<NonuniformTableLinear<double>> SatFuncGwseg;
std::vector<SatFuncGwseg> satfunc_;
std::vector<int> cell_to_func_; // = SATNUM - 1 std::vector<int> cell_to_func_; // = SATNUM - 1
std::vector<int> cell_to_func_imb_; std::vector<int> cell_to_func_imb_;
@@ -143,9 +148,6 @@ namespace Opm
std::vector<EPSTransforms> eps_transf_hyst_; std::vector<EPSTransforms> eps_transf_hyst_;
std::vector<SatHyst> sat_hyst_; std::vector<SatHyst> sat_hyst_;
typedef SatFuncSet Funcs;
inline const Funcs& funcForCell(const int cell) const;
template<class T> template<class T>
void initEPS(Opm::DeckConstPtr deck, void initEPS(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState, Opm::EclipseStateConstPtr eclipseState,
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
+110 -95
View File
@@ -24,6 +24,7 @@
#include <opm/core/utility/UniformTableLinear.hpp> #include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp> #include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp> #include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/simulator/ExplicitArraysFluidState.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/grid/GridHelpers.hpp> #include <opm/core/grid/GridHelpers.hpp>
@@ -123,12 +124,12 @@ namespace Opm
} }
} }
// Initialize tables. // Initialize saturation function objects.
satfuncset_.resize(num_tables); satfunc_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) { for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(eclipseState, table, phase_usage_, -1); satfunc_[table].init(eclipseState, table, phase_usage_, -1);
} }
// Check EHYSTR status // Check EHYSTR status
do_hyst_ = false; do_hyst_ = false;
if (hysteresis_switch && deck->hasKeyword("EHYSTR")) { if (hysteresis_switch && deck->hasKeyword("EHYSTR")) {
@@ -274,27 +275,31 @@ namespace Opm
{ {
assert(cells != 0); assert(cells != 0);
ExplicitArraysFluidState fluidState;
fluidState.setSaturationArray(s);
const int np = phase_usage_.num_phases; const int np = phase_usage_.num_phases;
if (dkrds) { if (dkrds) {
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
fluidState.setIndex(i);
if (do_hyst_) { if (do_hyst_) {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]])); satfunc_[cell_to_func_[cells[i]]].evalKrDeriv(fluidState, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]]));
} else if (do_eps_) { } else if (do_eps_) {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]])); satfunc_[cell_to_func_[cells[i]]].evalKrDeriv(fluidState, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]]));
} else { } else {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i); satfunc_[cell_to_func_[cells[i]]].evalKrDeriv(fluidState, kr + np*i, dkrds + np*np*i);
} }
} }
} else { } else {
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
if (do_hyst_) { if (do_hyst_) {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]])); satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]]));
} else if (do_eps_) { } else if (do_eps_) {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i, &(eps_transf_[cells[i]])); satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i, &(eps_transf_[cells[i]]));
} else { } else {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i); satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i);
} }
} }
} }
@@ -322,31 +327,34 @@ namespace Opm
{ {
assert(cells != 0); assert(cells != 0);
ExplicitArraysFluidState fluidState;
fluidState.setSaturationArray(s);
const int np = phase_usage_.num_phases; const int np = phase_usage_.num_phases;
if (dpcds) { if (dpcds) {
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
fluidState.setIndex(i);
if (do_eps_) { if (do_eps_) {
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i, &(eps_transf_[cells[i]])); satfunc_[cell_to_func_[cells[i]]].evalPcDeriv(fluidState, pc + np*i, dpcds + np*np*i, &(eps_transf_[cells[i]]));
} else { } else {
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i); satfunc_[cell_to_func_[cells[i]]].evalPcDeriv(fluidState, pc + np*i, dpcds + np*np*i);
} }
} }
} else { } else {
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
fluidState.setIndex(i);
if (do_eps_) { if (do_eps_) {
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i, &(eps_transf_[cells[i]])); satfunc_[cell_to_func_[cells[i]]].evalPc(fluidState, pc + np*i, &(eps_transf_[cells[i]]));
} else { } else {
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i); satfunc_[cell_to_func_[cells[i]]].evalPc(fluidState, pc + np*i);
} }
} }
} }
} }
/// Obtain the range of allowable saturation values. /// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices. /// \param[in] cells Array of n cell indices.
@@ -357,6 +365,19 @@ namespace Opm
const int* cells, const int* cells,
double* smin, double* smin,
double* smax) const double* smax) const
{
for (int cellIdx = 0; cellIdx < n; ++cellIdx) {
const SatFuncGwseg& satFunc = satfunc_[cell_to_func_[cellIdx]];
satRange_(satFunc, cellIdx, cells, smin, smax);
}
}
template <class SaturationFunction>
void SaturationPropsFromDeck::satRange_(const SaturationFunction& satFunc,
const int cellIdx,
const int* cells,
double* smin,
double* smax) const
{ {
assert(cells != 0); assert(cells != 0);
const int np = phase_usage_.num_phases; const int np = phase_usage_.num_phases;
@@ -365,35 +386,32 @@ namespace Opm
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid]; const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid];
const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
for (int i = 0; i < n; ++i) {
smin[np*i + opos] = 1.0; smin[np*cellIdx + opos] = 1.0;
smax[np*i + opos] = 1.0; smax[np*cellIdx + opos] = 1.0;
if (phase_usage_.phase_used[Aqua]) { if (phase_usage_.phase_used[BlackoilPhases::Aqua]) {
smin[np*i + wpos] = eps_transf_[cells[i]].wat.doNotScale ? funcForCell(cells[i]).smin_[wpos] smin[np*cellIdx + wpos] = eps_transf_[cells[cellIdx]].wat.doNotScale ? satFunc.smin_[wpos]
: eps_transf_[cells[i]].wat.smin; : eps_transf_[cells[cellIdx]].wat.smin;
smax[np*i + wpos] = eps_transf_[cells[i]].wat.doNotScale ? funcForCell(cells[i]).smax_[wpos] smax[np*cellIdx + wpos] = eps_transf_[cells[cellIdx]].wat.doNotScale ? satFunc.smax_[wpos]
: eps_transf_[cells[i]].wat.smax; : eps_transf_[cells[cellIdx]].wat.smax;
smin[np*i + opos] -= smax[np*i + wpos]; smin[np*cellIdx + opos] -= smax[np*cellIdx + wpos];
smax[np*i + opos] -= smin[np*i + wpos]; smax[np*cellIdx + opos] -= smin[np*cellIdx + wpos];
} }
if (phase_usage_.phase_used[Vapour]) { if (phase_usage_.phase_used[BlackoilPhases::Vapour]) {
smin[np*i + gpos] = eps_transf_[cells[i]].gas.doNotScale ? funcForCell(cells[i]).smin_[gpos] smin[np*cellIdx + gpos] = eps_transf_[cells[cellIdx]].gas.doNotScale ? satFunc.smin_[gpos]
: eps_transf_[cells[i]].gas.smin; : eps_transf_[cells[cellIdx]].gas.smin;
smax[np*i + gpos] = eps_transf_[cells[i]].gas.doNotScale ? funcForCell(cells[i]).smax_[gpos] smax[np*cellIdx + gpos] = eps_transf_[cells[cellIdx]].gas.doNotScale ? satFunc.smax_[gpos]
: eps_transf_[cells[i]].gas.smax; : eps_transf_[cells[cellIdx]].gas.smax;
smin[np*i + opos] -= smax[np*i + gpos]; smin[np*cellIdx + opos] -= smax[np*cellIdx + gpos];
smax[np*i + opos] -= smin[np*i + gpos]; smax[np*cellIdx + opos] -= smin[np*cellIdx + gpos];
} }
if (phase_usage_.phase_used[Vapour] && phase_usage_.phase_used[Aqua]) { if (phase_usage_.phase_used[BlackoilPhases::Vapour] && phase_usage_.phase_used[BlackoilPhases::Aqua]) {
smin[np*i + opos] = std::max(0.0,smin[np*i + opos]); smin[np*cellIdx + opos] = std::max(0.0,smin[np*cellIdx + opos]);
}
} }
} else { } else {
for (int i = 0; i < n; ++i) { for (int p = 0; p < np; ++p) {
for (int p = 0; p < np; ++p) { smin[np*cellIdx + p] = satFunc.smin_[p];
smin[np*i + p] = funcForCell(cells[i]).smin_[p]; smax[np*cellIdx + p] = satFunc.smax_[p];
smax[np*i + p] = funcForCell(cells[i]).smax_[p];
}
} }
} }
} }
@@ -413,7 +431,7 @@ namespace Opm
if (do_hyst_) { if (do_hyst_) {
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).updateSatHyst(s + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]])); satfunc_[cell_to_func_[cells[i]]].updateSatHyst(s + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]]));
} }
} }
} }
@@ -440,8 +458,11 @@ namespace Opm
const int max_np = BlackoilPhases::MaxNumPhases; const int max_np = BlackoilPhases::MaxNumPhases;
double s[max_np] = { 0.0 }; double s[max_np] = { 0.0 };
s[wpos] = swat; s[wpos] = swat;
ExplicitArraysFluidState fluidState;
fluidState.setSaturationArray(s);
fluidState.setIndex(0);
double pc[max_np] = { 0.0 }; double pc[max_np] = { 0.0 };
funcForCell(cell).evalPc(s, pc, &(eps_transf_[cell])); satfunc_[cell_to_func_[cell]].evalPc(fluidState, pc, &(eps_transf_[cell]));
if (pc[wpos] > pc_low_threshold) { if (pc[wpos] > pc_low_threshold) {
eps_transf_[cell].wat.pcFactor *= pcow/pc[wpos]; eps_transf_[cell].wat.pcFactor *= pcow/pc[wpos];
} }
@@ -452,13 +473,6 @@ namespace Opm
} }
// Map the cell number to the correct function set.
inline const SaturationPropsFromDeck::SatFuncSet&
SaturationPropsFromDeck::funcForCell(const int cell) const
{
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
// Initialize saturation scaling parameters // Initialize saturation scaling parameters
template<class T> template<class T>
void SaturationPropsFromDeck::initEPS(Opm::DeckConstPtr deck, void SaturationPropsFromDeck::initEPS(Opm::DeckConstPtr deck,
@@ -480,56 +494,57 @@ namespace Opm
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour]; const bool oilWater = phase_usage_.phase_used[BlackoilPhases::Aqua] && phase_usage_.phase_used[BlackoilPhases::Liquid] && !phase_usage_.phase_used[BlackoilPhases::Vapour];
const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; const bool oilGas = !phase_usage_.phase_used[BlackoilPhases::Aqua] && phase_usage_.phase_used[BlackoilPhases::Liquid] && phase_usage_.phase_used[BlackoilPhases::Vapour];
const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; const bool threephase = phase_usage_.phase_used[BlackoilPhases::Aqua] && phase_usage_.phase_used[BlackoilPhases::Liquid] && phase_usage_.phase_used[BlackoilPhases::Vapour];
for (int cell = 0; cell < number_of_cells; ++cell) { for (int cell = 0; cell < number_of_cells; ++cell) {
auto& satFunc = satfunc_[cell_to_func_[cell]];
if (threephase || oilWater) { if (threephase || oilWater) {
// ### krw // ### krw
initEPSParam(cell, eps_transf[cell].wat, false, initEPSParam(cell, eps_transf[cell].wat, false,
funcForCell(cell).smin_[wpos], satFunc.smin_[wpos],
funcForCell(cell).swcr_, satFunc.swcr_,
funcForCell(cell).smax_[wpos], satFunc.smax_[wpos],
funcForCell(cell).sowcr_, satFunc.sowcr_,
oilWater ? -1.0 : funcForCell(cell).smin_[gpos], oilWater ? -1.0 : satFunc.smin_[gpos],
funcForCell(cell).krwr_, satFunc.krwr_,
funcForCell(cell).krwmax_, satFunc.krwmax_,
funcForCell(cell).pcwmax_, satFunc.pcwmax_,
eps_vec[0], eps_vec[2], eps_vec[1], eps_vec[6], eps_vec[3], eps_vec[11], eps_vec[8], eps_vec[15]); eps_vec[0], eps_vec[2], eps_vec[1], eps_vec[6], eps_vec[3], eps_vec[11], eps_vec[8], eps_vec[15]);
// ### krow // ### krow
initEPSParam(cell, eps_transf[cell].watoil, true, initEPSParam(cell, eps_transf[cell].watoil, true,
0.0, 0.0,
funcForCell(cell).sowcr_, satFunc.sowcr_,
funcForCell(cell).smin_[wpos], satFunc.smin_[wpos],
funcForCell(cell).swcr_, satFunc.swcr_,
oilWater ? -1.0 : funcForCell(cell).smin_[gpos], oilWater ? -1.0 : satFunc.smin_[gpos],
funcForCell(cell).krorw_, satFunc.krorw_,
funcForCell(cell).kromax_, satFunc.kromax_,
0.0, 0.0,
eps_vec[0], eps_vec[6], eps_vec[0], eps_vec[2], eps_vec[3], eps_vec[13], eps_vec[10], dummy); eps_vec[0], eps_vec[6], eps_vec[0], eps_vec[2], eps_vec[3], eps_vec[13], eps_vec[10], dummy);
} }
if (threephase || oilGas) { if (threephase || oilGas) {
// ### krg // ### krg
initEPSParam(cell, eps_transf[cell].gas, false, initEPSParam(cell, eps_transf[cell].gas, false,
funcForCell(cell).smin_[gpos], satFunc.smin_[gpos],
funcForCell(cell).sgcr_, satFunc.sgcr_,
funcForCell(cell).smax_[gpos], satFunc.smax_[gpos],
funcForCell(cell).sogcr_, satFunc.sogcr_,
oilGas ? -1.0 : funcForCell(cell).smin_[wpos], oilGas ? -1.0 : satFunc.smin_[wpos],
funcForCell(cell).krgr_, satFunc.krgr_,
funcForCell(cell).krgmax_, satFunc.krgmax_,
funcForCell(cell).pcgmax_, satFunc.pcgmax_,
eps_vec[3], eps_vec[5], eps_vec[4], eps_vec[7], eps_vec[0], eps_vec[12], eps_vec[9], eps_vec[16]); eps_vec[3], eps_vec[5], eps_vec[4], eps_vec[7], eps_vec[0], eps_vec[12], eps_vec[9], eps_vec[16]);
// ### krog // ### krog
initEPSParam(cell, eps_transf[cell].gasoil, true, initEPSParam(cell, eps_transf[cell].gasoil, true,
0.0, 0.0,
funcForCell(cell).sogcr_, satFunc.sogcr_,
funcForCell(cell).smin_[gpos], satFunc.smin_[gpos],
funcForCell(cell).sgcr_, satFunc.sgcr_,
oilGas ? -1.0 : funcForCell(cell).smin_[wpos], oilGas ? -1.0 : satFunc.smin_[wpos],
funcForCell(cell).krorg_, satFunc.krorg_,
funcForCell(cell).kromax_, satFunc.kromax_,
0.0, 0.0,
eps_vec[3], eps_vec[7], eps_vec[3], eps_vec[5], eps_vec[0], eps_vec[14], eps_vec[10], dummy); eps_vec[3], eps_vec[7], eps_vec[3], eps_vec[5], eps_vec[0], eps_vec[14], eps_vec[10], dummy);
} }
@@ -547,9 +562,9 @@ namespace Opm
const std::string& keyword, const std::string& keyword,
std::vector<double>& scaleparam) std::vector<double>& scaleparam)
{ {
const bool useAqua = phase_usage_.phase_used[Aqua]; const bool useAqua = phase_usage_.phase_used[BlackoilPhases::Aqua];
const bool useLiquid = phase_usage_.phase_used[Liquid]; const bool useLiquid = phase_usage_.phase_used[BlackoilPhases::Liquid];
const bool useVapour = phase_usage_.phase_used[Vapour]; const bool useVapour = phase_usage_.phase_used[BlackoilPhases::Vapour];
bool useKeyword = deck->hasKeyword(keyword); bool useKeyword = deck->hasKeyword(keyword);
bool useStateKeyword = eclipseState->hasDoubleGridProperty(keyword); bool useStateKeyword = eclipseState->hasDoubleGridProperty(keyword);
const std::map<std::string, int> kw2tab = { const std::map<std::string, int> kw2tab = {
@@ -588,49 +603,49 @@ namespace Opm
itab = 1; itab = 1;
scaleparam.resize(number_of_cells); scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i) for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwmax_; scaleparam[i] = satfunc_[cell_to_func_[i]].krwmax_;
} }
} else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) { } else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 1))) { if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 1))) {
itab = 2; itab = 2;
scaleparam.resize(number_of_cells); scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i) for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgmax_; scaleparam[i] = satfunc_[cell_to_func_[i]].krgmax_;
} }
} else if (keyword == std::string("KRO") || keyword == std::string("IKRO") ) { } else if (keyword == std::string("KRO") || keyword == std::string("IKRO") ) {
if (useLiquid && (useKeyword || columnIsMasked_(deck, "ENKRVD", 2))) { if (useLiquid && (useKeyword || columnIsMasked_(deck, "ENKRVD", 2))) {
itab = 3; itab = 3;
scaleparam.resize(number_of_cells); scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i) for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).kromax_; scaleparam[i] = satfunc_[cell_to_func_[i]].kromax_;
} }
} else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) { } else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENKRVD", 3))) { if (useAqua && (useKeyword || columnIsMasked_(deck, "ENKRVD", 3))) {
itab = 4; itab = 4;
scaleparam.resize(number_of_cells); scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i) for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwr_; scaleparam[i] = satfunc_[cell_to_func_[i]].krwr_;
} }
} else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) { } else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 4))) { if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 4))) {
itab = 5; itab = 5;
scaleparam.resize(number_of_cells); scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i) for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgr_; scaleparam[i] = satfunc_[cell_to_func_[i]].krgr_;
} }
} else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) { } else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENKRVD", 5))) { if (useAqua && (useKeyword || columnIsMasked_(deck, "ENKRVD", 5))) {
itab = 6; itab = 6;
scaleparam.resize(number_of_cells); scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i) for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorw_; scaleparam[i] = satfunc_[cell_to_func_[i]].krorw_;
} }
} else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) { } else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 6))) { if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 6))) {
itab = 7; itab = 7;
scaleparam.resize(number_of_cells); scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i) for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorg_; scaleparam[i] = satfunc_[cell_to_func_[i]].krorg_;
} }
} else { } else {
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'"); OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
@@ -651,11 +666,11 @@ namespace Opm
if (useAqua && (keyword == std::string("PCW") || keyword == std::string("IPCW")) ) { if (useAqua && (keyword == std::string("PCW") || keyword == std::string("IPCW")) ) {
scaleparam.resize(number_of_cells); scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i) for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).pcwmax_; scaleparam[i] = satfunc_[cell_to_func_[i]].pcwmax_;
} else if (useVapour && (keyword == std::string("PCG") || keyword == std::string("IPCG")) ) { } else if (useVapour && (keyword == std::string("PCG") || keyword == std::string("IPCG")) ) {
scaleparam.resize(number_of_cells); scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i) for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).pcgmax_; scaleparam[i] = satfunc_[cell_to_func_[i]].pcgmax_;
} }
} }