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

remove the now obsolete SatFunc* classes

Browse Source
This commit is contained in:
Andreas Lauser committed 2015-09-02 12:29:18 +02:00
1 parent 921f377a1d
commit a1eb108528
8 files changed
-1972
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

CMakeLists_files.cmake
-6
View File
@@ -98,7 +98,6 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/props/rock/RockBasic.cpp
opm/core/props/rock/RockCompressibility.cpp
opm/core/props/rock/RockFromDeck.cpp
opm/core/props/satfunc/SatFuncSimple.cpp
opm/core/props/satfunc/SaturationPropsBasic.cpp
opm/core/simulator/AdaptiveSimulatorTimer.cpp
opm/core/simulator/BlackoilState.cpp
@@ -363,11 +362,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/props/rock/RockBasic.hpp
opm/core/props/rock/RockCompressibility.hpp
opm/core/props/rock/RockFromDeck.hpp
opm/core/props/satfunc/SatFuncBase.hpp
opm/core/props/satfunc/SatFuncGwseg.hpp
opm/core/props/satfunc/SatFuncMultiplexer.hpp
opm/core/props/satfunc/SatFuncSimple.hpp
opm/core/props/satfunc/SatFuncStone2.hpp
opm/core/props/satfunc/SaturationPropsBasic.hpp
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
opm/core/props/satfunc/SatFuncBase.hpp
-329
View File
@@ -1,329 +0,0 @@
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef SATFUNCBASE_HPP
#define SATFUNCBASE_HPP
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
namespace Opm
{
// Transforms for saturation table scaling
struct EPSTransforms {
struct Transform {
bool doNotScale;
bool do_3pt;
double smin;
double scr;
double sr;
double smax;
double slope1;
double slope2;
double scaleSat(double ss, double s_r, double s_cr, double s_max) const;
double scaleSatInv(double s, double s_r, double s_cr, double s_max) const;
double scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const; // Returns scaleSat'(s)
double scaleSatPc(double s, double s_min, double s_max) const;
double scaleSatDerivPc(double s, double s_min, double s_max) const; // Returns scaleSatPc'(s)
bool doKrMax;
bool doKrCrit;
bool doSatInterp;
double krsr;
double krmax;
double krSlopeMax;
double krSlopeCrit;
double scaleKr(double s, double kr, double krsr_tab) const;
double scaleKrDeriv(double s, double krDeriv) const; // Returns scaleKr'(kr(scaleSat(s)))*kr'((scaleSat(s))
double pcFactor; // Scaling factor for capillary pressure.
void printMe(std::ostream & out);
};
Transform wat;
Transform watoil;
Transform gas;
Transform gasoil;
};
// Hysteresis
struct SatHyst {
SatHyst() {
sg_hyst = -1.0;
sow_hyst = -1.0;
sg_shift = 0.0;
sow_shift = 0.0;
}
double sg_hyst;
double sg_shift;
double sow_hyst;
double sow_shift;
void printMe(std::ostream & out);
};
template <class TableType>
class SatFuncBase : public BlackoilPhases
{
public:
virtual ~SatFuncBase() {};
void init(Opm::EclipseStateConstPtr eclipseState,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void updateSatHyst(const double* s,
const EPSTransforms* epst,
const EPSTransforms* epst_hyst,
SatHyst* sat_hyst) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double krgmax_; // Max gas relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double sgcr_; // Critical gas saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double krgr_; // Gas relperm at critical oil-in-gas saturation.
double sowcr_; // Critical oil-in-water saturation.
double sogcr_; // Critical oil-in-gas-and-connate-water saturation.
double krorw_; // Oil relperm at critical water saturation.
double krorg_; // Oil relperm at critical gas saturation.
double pcwmax_; // Max oil-water capillary pressure.
double pcgmax_; // Max gas-oil capillary pressure.
protected:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
TableType krw_;
TableType krow_;
TableType pcow_;
TableType krg_;
TableType krog_; TableType pcog_;
double krocw_; // = krow_(s_wc)
private:
void extendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
void initializeTableType(TableType& table,
const std::vector<double>& arg,
const std::vector<double>& value,
const int samples);
};
template <class TableType>
void SatFuncBase<TableType>::init(Opm::EclipseStateConstPtr eclipseState,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
{
std::shared_ptr<const TableManager> tables = eclipseState->getTableManager();
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const Opm::SwofTable& swof(tables->getSwofTables()[table_num]);
const std::vector<double>& sw = swof.getSwColumn();
const std::vector<double>& krw = swof.getKrwColumn();
const std::vector<double>& krow = swof.getKrowColumn();
const std::vector<double>& pcow = swof.getPcowColumn();
if (krw.front() != 0.0 || krow.back() != 0.0) {
OPM_THROW(std::runtime_error, "Error SWOF data - non-zero krw(swco) and/or krow(1-sor)");
}
// Extend the tables with constant values such that the
// derivatives at the endpoints are zero
int n = sw.size();
std::vector<double> sw_ex(n+2);
std::vector<double> krw_ex(n+2);
std::vector<double> krow_ex(n+2);
std::vector<double> pcow_ex(n+2);
extendTable(sw,sw_ex,1);
extendTable(krw,krw_ex,0);
extendTable(krow,krow_ex,0);
extendTable(pcow,pcow_ex,0);
initializeTableType(krw_,sw_ex, krw_ex, samples);
initializeTableType(krow_,sw_ex, krow_ex, samples);
initializeTableType(pcow_,sw_ex, pcow_ex, samples);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
krwmax_ = krw.back();
kromax_ = krow.front();
swcr_ = swmax;
sowcr_ = 1.0 - swco;
krwr_ = krw.back();
krorw_ = krow.front();
for (std::vector<double>::size_type i=1; i<sw.size(); ++i) {
if (krw[i]> 0.0) {
swcr_ = sw[i-1];
krorw_ = krow[i-1];
break;
}
}
for (std::vector<double>::size_type i=sw.size()-1; i>=1; --i) {
if (krow[i-1]> 0.0) {
sowcr_ = 1.0 - sw[i];
krwr_ = krw[i];
break;
}
}
pcwmax_ = pcow.front();
}
if (phase_usage.phase_used[Vapour]) {
const Opm::SgofTable& sgof = tables->getSgofTables()[table_num];
const std::vector<double>& sg = sgof.getSgColumn();
const std::vector<double>& krg = sgof.getKrgColumn();
const std::vector<double>& krog = sgof.getKrogColumn();
const std::vector<double>& pcog = sgof.getPcogColumn();
// Extend the tables with constant values such that the
// derivatives at the endpoints are zero
int n = sg.size();
std::vector<double> sg_ex(n+2);
std::vector<double> krg_ex(n+2);
std::vector<double> krog_ex(n+2);
std::vector<double> pcog_ex(n+2);
extendTable(sg,sg_ex,1);
extendTable(krg,krg_ex,0);
extendTable(krog,krog_ex,0);
extendTable(pcog,pcog_ex,0);
initializeTableType(krg_,sg_ex, krg_ex, samples);
initializeTableType(krog_,sg_ex, krog_ex, samples);
initializeTableType(pcog_,sg_ex, pcog_ex, samples);
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
smin_[phase_usage.phase_pos[Vapour]] = sg.front();
krgmax_ = krg.back();
sgcr_ = sg.front();
sogcr_ = 1.0 - sg.back();
krgr_ = krg.back();
krorg_ = krg.front();
for (std::vector<double>::size_type i=1; i<sg.size(); ++i) {
if (krg[i]> 0.0) {
sgcr_ = sg[i-1];
krorg_ = krog[i-1];
break;
}
}
for (std::vector<double>::size_type i=sg.size()-1; i>=1; --i) {
if (krog[i-1]> 0.0) {
sogcr_ = 1.0 - sg[i];
krgr_ = krg[i];
break;
}
}
pcgmax_ = pcog.back();
}
if (phase_usage.phase_used[Vapour] && phase_usage.phase_used[Aqua]) {
sowcr_ -= smin_[phase_usage.phase_pos[Vapour]];
sogcr_ -= smin_[phase_usage.phase_pos[Aqua]];
smin_[phase_usage.phase_pos[Liquid]] = 0.0;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]]
- smin_[phase_usage.phase_pos[Vapour]]; // First entry in SGOF-table supposed to be zero anyway ...
} else if (phase_usage.phase_used[Aqua]) {
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Aqua]];
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]];
} else if (phase_usage.phase_used[Vapour]) {
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Vapour]];
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Vapour]];
}
}
template <class TableType>
void SatFuncBase<TableType>::updateSatHyst(const double* s,
const EPSTransforms* epst,
const EPSTransforms* epst_hyst,
SatHyst* sat_hyst) const
{
if (phase_usage.phase_used[Aqua] && phase_usage.phase_used[Vapour]) { //Water/Oil/Gas
int opos = phase_usage.phase_pos[BlackoilPhases::Liquid];
int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour];
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
if (s[opos] > sat_hyst->sow_hyst)
{
sat_hyst->sow_hyst = s[opos];
double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]);
double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]);
sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst;
}
if (s[gpos] > sat_hyst->sg_hyst)
{
sat_hyst->sg_hyst = s[gpos];
double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]);
double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]);
sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst;
}
} else if (phase_usage.phase_used[Aqua]) { //Water/oil
int opos = phase_usage.phase_pos[BlackoilPhases::Liquid];
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
if (s[opos] > sat_hyst->sow_hyst)
{
sat_hyst->sow_hyst = s[opos];
double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]);
double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]);
sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst;
}
} else if (phase_usage.phase_used[Vapour]) {//Gas/Oil
int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour];
if (s[gpos] > sat_hyst->sg_hyst)
{
sat_hyst->sg_hyst = s[gpos];
double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]);
double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]);
sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst;
}
}
}
template <class TableType>
void SatFuncBase<TableType>::extendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const
{
int n = xv.size();
xv_ex[0] = xv[0]-pm;
xv_ex[n+1] = xv[n-1]+pm;
for (int i=0; i<n; i++)
{
xv_ex[i+1] = xv[i];
}
}
} // namespace Opm
#endif // SATFUNCBASE_HPP
opm/core/props/satfunc/SatFuncGwseg.hpp
-637
View File
@@ -1,637 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef SATFUNCGWSEG_HPP
#define SATFUNCGWSEG_HPP
#include <opm/core/props/satfunc/SatFuncBase.hpp>
namespace Opm
{
template<class TableType>
class SatFuncGwseg : public SatFuncBase<TableType>
{
public:
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr) const;
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const;
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc) const;
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const;
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<NonuniformTableLinear<double> > SatFuncGwsegNonuniform;
template<class TableType>
template <class FluidState>
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr) const
{
if (this->phase_usage.num_phases == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-5;
swco = std::min(swco,sw-eps);
// xw and xg are the fractions occupied by water and gas zones.
const double ssg = sw - swco + sg;
const double xw = (sw - swco) / ssg;
const double xg = 1 - xw;
const double ssw = sg + sw;
const double krw = this->krw_(sw);
const double krg = this->krg_(sg);
const double krow = this->krow_(ssw);
const double krog = this->krog_(ssg);
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
template <class FluidState>
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst) const
{
if (this->phase_usage.num_phases == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
// TODO Also consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin;
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
const double ssg = std::max(sg + sw - swco, eps);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
double ssow = 1.0-ssw;
double ssog = 1.0-ssg-swco;
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
const double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
const double krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
const double krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
const double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
// xw and xg are the fractions occupied by water and gas zones.
const double xw = (sw - swco) / d;
const double xg = 1 - xw;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
return;
}
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
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) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
// TODO Consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin;
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
const double ssg = std::max(sg + sw - swco, eps);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
double ssow = 1.0-ssw;
double ssog = 1.0-ssg-swco;
// The code below corresponds to EHYSTR * 0 * * KR/
// - wettability properties water>oil>gas.
// - Carlsen hysteresis model for non-wetting (scanning=shifted_imb). No hysteresis for wetting phase.
// The imb-curve currently only differs from drainage curves via endpoint scaling ...
// Water - use drainage curve only
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
// Gas
double krg;
if (sg >= sat_hyst->sg_hyst) { // Drainage
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
} else { // Imbibition
double sg_shifted = sg + sat_hyst->sg_shift;
double _sg = epst_hyst->gas.scaleSat(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
krg = epst_hyst->gas.scaleKr(sg_shifted, this->krg_(_sg), this->krgr_);
}
// Oil in water
double krow;
if (ssow >= sat_hyst->sow_hyst) { // Drainage
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
} else { // Imbibition
double ssow_shifted = ssow + sat_hyst->sow_shift;
double _ssow = epst_hyst->watoil.scaleSat(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(1.0-_ssow), this->krorw_);
}
// Oil in gas and connate water - use drainage curve only
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
// xw and xg are the fractions occupied by water and gas zones.
const double xw = (sw - swco) / d;
const double xg = 1 - xw;
// relperms
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
return;
}
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
template<class FluidState>
void SatFuncGwseg<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-5;
swco = std::min(swco,sw-eps);
// xw and xg are the fractions occupied by water and gas zones.
const double ssw = sg + sw;
// d = ssg = sw - swco + sg (using 'd' for consistency with mrst docs).
const double d = sg + sw - swco;
const double xw = (sw - swco) / d;
const double krw = this->krw_(sw);
const double krg = this->krg_(sg);
const double krow = this->krow_(ssw);
const double krog = this->krog_(d);
const double xg = 1 - xw;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
// Derivatives.
const double dkrww = this->krw_.derivative(sw);
const double dkrgg = this->krg_.derivative(sg);
const double dkrow = this->krow_.derivative(ssw);
const double dkrog = this->krog_.derivative(d);
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
double dkrow = this->krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
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;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
// TODO Also consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin;
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
const double ssg = std::max(sg + sw - swco, eps);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
double ssow = 1.0-ssw;
double ssog = 1.0-ssg-swco;
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _dsdsw = epst->wat.scaleSatDeriv(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _dsdsg = epst->gas.scaleSatDeriv(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssow = epst->watoil.scaleSatDeriv(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssog = epst->gasoil.scaleSatDeriv(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
const double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
const double krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
const double krow = epst->watoil.scaleKr(ssow, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
const double krog = epst->gasoil.scaleKr(ssog, this->krog_(std::max(1.0-_ssog-_swco,eps)), this->krorg_);
// xw and xg are the fractions occupied by water and gas zones.
const double xw = (sw - swco) / d;
const double xg = 1 - xw;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
// Derivatives.
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(sw, this->krw_.derivative(_sw));
double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(sg, this->krg_.derivative(_sg));
double dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(std::max(1.0-_ssog-_swco,eps)));
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
return;
}
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
double dkrow = this->krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
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;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
// TODO Also consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin;
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
const double ssg = std::max(sg + sw - swco, eps);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
double ssow = 1.0-ssw;
double ssog = 1.0-ssg-swco;
// The code below corresponds to EHYSTR * 0 * * KR/
// - wettability properties water>oil>gas.
// - Carlsen hysteresis model for non-wetting (scanning=shifted_imb). No hysteresis for wetting phase.
// The imb-curve currently only differs from drainage curves via endpoint scaling ...
// Water - use drainage curve only
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _dsdsw = epst->wat.scaleSatDeriv(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(sw, this->krw_.derivative(_sw));
// Gas
double krg, dkrgg;
if (sg >= sat_hyst->sg_hyst) { // Drainage
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _dsdsg = epst->gas.scaleSatDeriv(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
dkrgg = _dsdsg*epst->gas.scaleKrDeriv(sg, this->krg_.derivative(_sg));
} else { // Imbibition
double sg_shifted = sg + sat_hyst->sg_shift;
double _sg = epst_hyst->gas.scaleSat(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _dsdsg = epst_hyst->gas.scaleSatDeriv(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
krg = epst_hyst->gas.scaleKr(sg_shifted, this->krg_(_sg), this->krgr_);
dkrgg = _dsdsg*epst_hyst->gas.scaleKrDeriv(sg_shifted, this->krg_.derivative(_sg));
}
// Oil in water
double krow, dkrow;
if (ssow >= sat_hyst->sow_hyst) { // Drainage
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssow = epst->watoil.scaleSatDeriv(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
krow = epst->watoil.scaleKr(ssow, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
} else { // Imbibition
double ssow_shifted = ssow + sat_hyst->sow_shift;
double _ssow = epst_hyst->watoil.scaleSat(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssow = epst_hyst->watoil.scaleSatDeriv(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
dkrow = _dsdssow*epst_hyst->watoil.scaleKrDeriv(ssow_shifted, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
}
// Oil in gas and connate water - use drainage curve only
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssog = epst->gasoil.scaleSatDeriv(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double krog = epst->gasoil.scaleKr(ssog, this->krog_(std::max(1.0-_ssog-_swco,eps)), this->krorg_);
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(std::max(1.0-_ssog-_swco,eps)));
// xw and xg are the fractions occupied by water and gas zones.
const double xw = (sw - swco) / d;
const double xg = 1 - xw;
// relperms
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
// Derivatives.
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
return;
}
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
double dkrow = this->krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
template <class FluidState>
void SatFuncGwseg<TableType>::evalPc(const FluidState& fluidState, double* pc) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(fluidState.saturation(pos));
}
}
template<class TableType>
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;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
double _sw = epst->wat.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->wat.pcFactor*this->pcow_(_sw);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
double _sg = epst->gas.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->gas.pcFactor*this->pcog_(_sg);
}
}
template<class TableType>
template <class FluidState>
void SatFuncGwseg<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = this->phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcow_.derivative(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcog_.derivative(fluidState.saturation(pos));
}
}
template<class TableType>
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
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = this->phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
double _sw = epst->wat.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->wat.pcFactor*this->pcow_(_sw);
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);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
double _sg = epst->gas.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->gas.pcFactor*this->pcog_(_sg);
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);
}
}
} // namespace Opm
#endif // SATFUNCGWSEG_HPP
opm/core/props/satfunc/SatFuncMultiplexer.hpp
-268
View File
@@ -1,268 +0,0 @@
/*
Copyright 2015 Andreas Lauser
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
/*!
* \file
* \copydoc Opm::SatFuncMultiplexer
*/
#ifndef OPM_SAT_FUNC_MULTIPLEXER_HPP
#define OPM_SAT_FUNC_MULTIPLEXER_HPP
#include "SatFuncGwseg.hpp"
#include "SatFuncSimple.hpp"
#include "SatFuncStone2.hpp"
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <memory>
namespace Opm
{
/*!
* \brief This is a multiplexer class which allows to saturation function to be used at
* runtime.
*/
class SatFuncMultiplexer
{
public:
typedef NonuniformTableLinear<double> TableType;
enum SatFuncType {
Gwseg, // <- supposed to be the default
Stone2,
Simple
};
// this is a helper macro which helps to save us from RSI in the following
#define OPM_MULTIPLEXER_CONST const
#define OPM_MULTIPLEXER_NON_CONST
#define OPM_MULTIPLEXER_SATFUNC_CALL(codeToCall, CONST) \
switch (satFuncType_) { \
case Gwseg: { \
typedef SatFuncGwseg<TableType> SatFunc; \
__attribute__((unused)) CONST SatFunc& satFunc = \
*static_cast<SatFunc*>(satFunc_.get()); \
codeToCall; \
break; \
} \
case Stone2: { \
typedef SatFuncStone2<TableType> SatFunc; \
__attribute__((unused)) CONST SatFunc& satFunc = \
*static_cast<SatFunc*>(satFunc_.get()); \
codeToCall; \
break; \
} \
case Simple: { \
typedef SatFuncSimple<TableType> SatFunc; \
__attribute__((unused)) CONST SatFunc& satFunc = \
*static_cast<SatFunc*>(satFunc_.get()); \
codeToCall; \
break; \
} \
};
SatFuncMultiplexer()
{}
~SatFuncMultiplexer()
{}
// this constructor not really a copy constructor, but it is
// required to make std::unique_ptr happy
SatFuncMultiplexer(const SatFuncMultiplexer&)
{}
// this operator does not do anything and is thus not a copy operator, but it is
// required to make std::unique_ptr happy on old compilers
SatFuncMultiplexer& operator=(const SatFuncMultiplexer&)
{ return *this; }
/*!
* \brief Pick the correct saturation function type and initialize the object using
* an ECL deck.
*/
void initFromDeck(Opm::EclipseStateConstPtr eclState,
size_t tableIdx,
SatFuncType satFuncType)
{
auto phaseUsage = phaseUsageFromDeck(eclState);
setSatFuncType(satFuncType);
OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.init(eclState,
tableIdx,
phaseUsage,
/*samples=*/0), OPM_MULTIPLEXER_NON_CONST);
}
/*!
* \brief Returns the type of the actually selected saturation function.
*/
void setSatFuncType(SatFuncType newSatFuncType)
{
// set the type of the saturation function
satFuncType_ = newSatFuncType;
// call the default constructor for the selected saturation function
switch (satFuncType_) {
case Gwseg: {
typedef SatFuncGwseg<TableType> SatFunc;
satFunc_.reset(new SatFunc);
break;
}
case Stone2: {
typedef SatFuncStone2<TableType> SatFunc;
satFunc_.reset(new SatFunc);
break;
}
case Simple: {
typedef SatFuncSimple<TableType> SatFunc;
satFunc_.reset(new SatFunc);
break;
}
}
}
/*!
* \brief Returns which saturation function is actually selected.
*/
SatFuncType satFuncType() const
{ return satFuncType_; }
/*!
* \brief Returns a pointer to the base saturation function.
*
* All underlying saturation functions derive from this class. It stores the stuff
* which is common for all saturation functions (e.g. the endpoint scaling
* parameters).
*/
const SatFuncBase<TableType>& getSatFuncBase() const
{ return *satFunc_; }
SatFuncBase<TableType>& getSatFuncBase()
{ return *satFunc_; }
/*!
* \brief Return the raw "Gwseg" saturation function object.
*
* If the saturation function type is not "Gwseg", an assertation is triggered for
* debug builds.
*/
SatFuncGwseg<TableType>& getGwseg()
{
assert(satFuncType_ == Gwseg);
return *static_cast<SatFuncGwseg<TableType>*>(satFunc_.get());
}
const SatFuncGwseg<TableType>& getGwseg() const
{
assert(satFuncType_ == Gwseg);
return *static_cast<const SatFuncGwseg<TableType>*>(satFunc_.get());
}
/*!
* \brief Return the raw "Simple" saturation function object.
*
* If the saturation function type is not "Simple", an assertation is triggered for
* debug builds.
*/
SatFuncSimple<TableType>& getSimple()
{
assert(satFuncType_ == Simple);
return *static_cast<SatFuncSimple<TableType>*>(satFunc_.get());
}
const SatFuncSimple<TableType>& getSimple() const
{
assert(satFuncType_ == Simple);
return *static_cast<const SatFuncSimple<TableType>*>(satFunc_.get());
}
/*!
* \brief Return the raw "Stone2" saturation function object.
*
* If the saturation function type is not "Stone2", an assertation is triggered for
* debug builds.
*/
SatFuncStone2<TableType>& getStone2()
{
assert(satFuncType_ == Stone2);
return *static_cast<SatFuncStone2<TableType>*>(satFunc_.get());
}
const SatFuncStone2<TableType>& getStone2() const
{
assert(satFuncType_ == Stone2);
return *static_cast<const SatFuncStone2<TableType>*>(satFunc_.get());
}
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKr(fluidState, kr), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKrDeriv(fluidState, kr, dkrds), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPc(fluidState, pc), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPcDeriv(fluidState, pc, dpcds), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKr(fluidState, kr, epst), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKr(fluidState, kr, epst, epst_hyst, sat_hyst), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKrDeriv(fluidState, kr, dkrds, epst), OPM_MULTIPLEXER_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
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKrDeriv(fluidState, kr, dkrds, epst, epst_hyst, sat_hyst), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc, const EPSTransforms* epst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPc(fluidState, pc, epst), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds, const EPSTransforms* epst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPcDeriv(fluidState, pc, dpcds, epst), OPM_MULTIPLEXER_CONST); }
// undefine the helper macro here because we don't need it anymore
#undef OPM_MULTIPLEXER_SATFUNC_CALL
#undef OPM_MULTIPLEXER_NON_CONST
#undef OPM_MULTIPLEXER_CONST
private:
SatFuncType satFuncType_;
std::unique_ptr<SatFuncBase<TableType> > satFunc_;
};
} // namespace Opm
#endif // OPM_SAT_FUNC_MULTIPLEXER_HPP
opm/core/props/satfunc/SatFuncSimple.cpp
-223
View File
@@ -1,223 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/props/satfunc/SatFuncSimple.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>
namespace Opm
{
// SatFuncSimple.cpp can now be removed and the code below moved to SatFuncBase.cpp ...
template<>
void SatFuncBase<NonuniformTableLinear<double> >::initializeTableType(NonuniformTableLinear<double> & table,
const std::vector<double>& arg,
const std::vector<double>& value,
const int /* samples */)
{
table = NonuniformTableLinear<double>(arg, value);
}
template<>
void SatFuncBase<UniformTableLinear<double> >::initializeTableType(UniformTableLinear<double> & table,
const std::vector<double>& arg,
const std::vector<double>& value,
const int samples)
{
buildUniformMonotoneTable(arg, value, samples, table);
}
double EPSTransforms::Transform::scaleSat(double s, double s_r, double s_cr, double s_max) const
{
if (doNotScale) {
return s;
} else if (!do_3pt) { // 2-pt
if (s <= scr) {
return s_cr;
} else {
return (s >= smax) ? s_max : s_cr + (s-scr)*slope1;
}
} else if (s <= sr) {
return (s <= scr) ? s_cr : s_cr+(s-scr)*slope1;
} else {
return (s >= smax) ? s_max : s_r+(s-sr)*slope2;
}
}
double EPSTransforms::Transform::scaleSatInv(double ss, double s_r, double s_cr, double s_max) const
{
if (doNotScale) {
return ss;
} else if (!do_3pt) { // 2-pt
if (ss <= s_cr) {
return scr;
} else {
return (ss >= s_max) ? smax : scr + (ss-s_cr)/slope1;
}
} else if (ss <= s_r) {
return (ss <= s_cr) ? scr : scr+(ss-s_cr)/slope1;
} else {
return (ss >= s_max) ? smax : sr+(ss-s_r)/slope2;
}
}
double EPSTransforms::Transform::scaleSatDeriv(double s, double /* s_r */, double /* s_cr */, double /* s_max */) const
{
if (doNotScale) {
return 1.0;
} else if (!do_3pt) { // 2-pt
if (s <= scr) {
return 0.0;
} else {
return (s >= smax) ? 0.0 : slope1;
}
} else if (s <= sr) {
return (s <= scr) ? 0.0 : slope1;
} else {
return (s >= smax) ? 0.0 : slope2;
}
}
double EPSTransforms::Transform::scaleSatPc(double s, double s_min, double s_max) const
{
if (doNotScale) {
return s;
} else if (s<=smin) {
return s_min;
} else if (s <= smax) {
return s_min + (s-smin)*(s_max-s_min)/(smax-smin);
} else {
return s_max;
}
}
double EPSTransforms::Transform::scaleSatDerivPc(double s, double s_min, double s_max) const
{
if (doNotScale) {
return 1.0;
} else if (s<smin) {
return 0.0;
} else if (s <= smax) {
return (s_max-s_min)/(smax-smin);
} else {
return 0.0;
}
}
double EPSTransforms::Transform::scaleKr(double s, double kr, double krsr_tab) const
{
if (doKrCrit) {
if (s <= scr) {
return 0.0;
} else if (s <= sr) {
return kr*krSlopeCrit;
} else if (s <= smax) {
if (doSatInterp)
return krsr + (s-sr)*krSlopeMax; // Note: Scaling independent of kr-value ...
else
return krsr + (kr-krsr_tab)*krSlopeMax;
} else {
return krmax;
}
} else if (doKrMax) {
if (s <= scr) {
return 0.0;
} else if (s <= smax) {
return kr*krSlopeMax;
} else {
return krmax;
}
} else {
return kr;
}
}
double EPSTransforms::Transform::scaleKrDeriv(double s, double krDeriv) const
{
if (doKrCrit) {
if (s < scr) {
return 0.0;
} else if (s <= sr) {
return krDeriv*krSlopeCrit;
} else if (s <= smax) {
if (doSatInterp)
return krSlopeMax; // Note: Scaling independent of kr-value ...
else
return krDeriv*krSlopeMax;
} else {
return 0.0;
}
} else if (doKrMax) {
if (s < scr) {
return 0.0;
} else if (s <= smax) {
return krDeriv*krSlopeMax;
} else {
return 0.0;
}
} else {
if (s < scr) {
return 0.0;
} else if (s <= smax) {
return krDeriv;
} else {
return 0.0;
}
}
}
void EPSTransforms::Transform::printMe(std::ostream & out)
{
out << "doNotScale: " << doNotScale << std::endl;
out << "do_3pt: " << do_3pt << std::endl;
out << "smin: " << smin << std::endl;
out << "scr: " << scr << std::endl;
out << "sr: " << sr << std::endl;
out << "smax: " << smax << std::endl;
out << "slope1: " << slope1 << std::endl;
out << "slope2: " << slope2 << std::endl;
out << "doKrMax: " << doKrMax << std::endl;
out << "doKrCrit: " << doKrCrit << std::endl;
out << "doSatInterp: " << doSatInterp << std::endl;
out << "krsr: " << krsr << std::endl;
out << "krmax: " << krmax << std::endl;
out << "krSlopeMax: " << krSlopeMax << std::endl;
out << "krSlopeCrit: " << krSlopeCrit << std::endl;
}
void SatHyst::printMe(std::ostream & out)
{
out << "sg_hyst: " << sg_hyst << std::endl;
out << "sg_shift: " << sg_shift << std::endl;
out << "sow_hyst: " << sow_hyst << std::endl;
out << "sow_shift: " << sow_shift << std::endl;
}
} // namespace Opm
opm/core/props/satfunc/SatFuncSimple.hpp
-297
View File
@@ -1,297 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef SATFUNCSIMPLE_HPP
#define SATFUNCSIMPLE_HPP
#include <opm/core/props/satfunc/SatFuncBase.hpp>
namespace Opm
{
template<class TableType>
class SatFuncSimple : public SatFuncBase<TableType>
{
public:
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr) const;
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const;
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc) const;
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const;
template <class FluidState>
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
template <class FluidState>
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
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
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
template <class FluidState>
void evalPc(const FluidState& /* fluidState */, double* /* pc */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
template <class FluidState>
void evalPcDeriv(const FluidState& /* fluidState */, double* /* pc */, double* /* dpcds */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
private:
};
typedef SatFuncSimple<UniformTableLinear<double> > SatFuncSimpleUniform;
typedef SatFuncSimple<NonuniformTableLinear<double> > SatFuncSimpleNonuniform;
template<class TableType>
template<class FluidState>
void SatFuncSimple<TableType>::evalKr(const FluidState& fluidState, double* kr) const
{
if (this->phase_usage.num_phases == 3) {
// A simplified relative permeability model.
double sw = fluidState.saturation(BlackoilPhases::Aqua);
double sg = fluidState.saturation(BlackoilPhases::Vapour);
double krw = this->krw_(sw);
double krg = this->krg_(sg);
double krow = this->krow_(sw + sg); // = 1 - so
// double krog = krog_(sg); // = 1 - so - sw
// double krocw = krocw_;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = krow;
if (kr[BlackoilPhases::Liquid] < 0.0) {
kr[BlackoilPhases::Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double so = fluidState.saturation(opos);
double krow = this->krow_(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
template<class FluidState>
void SatFuncSimple<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// A simplified relative permeability model.
double sw = fluidState.saturation(BlackoilPhases::Aqua);
double sg = fluidState.saturation(BlackoilPhases::Vapour);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krow = this->krow_(sw + sg);
double dkrow = this->krow_.derivative(sw + sg);
// double krog = krog_(sg);
// double dkrog = krog_.derivative(sg);
// double krocw = krocw_;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = krow;
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[BlackoilPhases::Liquid] < 0.0) {
kr[BlackoilPhases::Liquid] = 0.0;
}
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
//dkrds[Liquid + Aqua*np] = dkrow;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Liquid*np] = -dkrow;
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = 0.0;
//krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
//+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double so = fluidState.saturation(opos);
double krow = this->krow_(1.0-so);
double dkrow = this->krow_.derivative(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
template<class FluidState>
void SatFuncSimple<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
// A simplified relative permeability model.
// Define KR(s) = scaleKr(kr(scalSat(s)))
// Thus KR'(s) = scaleKr'(kr(scaleSat(s)))*kr'((scaleSat(s))*scaleSat'(s)
double _sw = epst->wat.scaleSat(fluidState.saturation(wpos), 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _dsdsw = epst->wat.scaleSatDeriv(fluidState.saturation(wpos), 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _sg = epst->gas.scaleSat(fluidState.saturation(gpos), 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _dsdsg = epst->gas.scaleSatDeriv(fluidState.saturation(gpos), 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _sow = epst->watoil.scaleSat(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdsow = epst->watoil.scaleSatDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
//double _sog = epst->gasoil.scaleSat(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
//double _dsdsog = epst->gasoil.scaleSatDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double krw = epst->wat.scaleKr(fluidState.saturation(wpos), this->krw_(_sw), this->krwr_);
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(fluidState.saturation(wpos), this->krw_.derivative(_sw));
double krg = epst->gas.scaleKr(fluidState.saturation(gpos), this->krg_(_sg), this->krgr_);
double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(fluidState.saturation(gpos), this->krg_.derivative(_sg));
// TODO Check the arguments to the krow- and krog-tables below...
double krow = epst->watoil.scaleKr(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krow_(1.0-_sow-this->smin_[gpos]), this->krorw_); // ????
double dkrow = _dsdsow*epst->watoil.scaleKrDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krow_.derivative(1.0-_sow-this->smin_[gpos])); // ????
//double krog = epst->gasoil.scaleKr(this->krog_(1.0-_sog-this->smin_[wpos]), 1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krorg_); // ????
//double dkrog = _dsdsog*epst->gasoil.scaleKrDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krog_.derivative(1.0-_sog-this->smin_[wpos])); // ????
// double krocw = krocw_;
kr[wpos] = krw;
kr[gpos] = krg;
kr[BlackoilPhases::Liquid] = krow;
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[BlackoilPhases::Liquid] < 0.0) {
kr[BlackoilPhases::Liquid] = 0.0;
}
dkrds[wpos + wpos*np] = dkrww;
dkrds[gpos + gpos*np] = dkrgg;
//dkrds[Liquid + Aqua*np] = dkrow;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Liquid*np] = -dkrow;
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[BlackoilPhases::Liquid + gpos*np] = 0.0;
//krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
//+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double so = fluidState.saturation(opos);
double krow = this->krow_(1.0-so);
double dkrow = this->krow_.derivative(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
template<class FluidState>
void SatFuncSimple<TableType>::evalPc(const FluidState& fluidState, double* pc) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(fluidState.saturation(pos));
}
}
template<class TableType>
template<class FluidState>
void SatFuncSimple<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = this->phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcow_.derivative(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcog_.derivative(fluidState.saturation(pos));
}
}
} // namespace Opm
#endif // SATFUNCSIMPLE_HPP
opm/core/props/satfunc/SatFuncStone2.hpp
-211
View File
@@ -1,211 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef SATFUNCSTONE2_HPP
#define SATFUNCSTONE2_HPP
#include <opm/core/props/satfunc/SatFuncBase.hpp>
namespace Opm
{
template<class TableType>
class SatFuncStone2 : public SatFuncBase<TableType>
{
public:
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr) const;
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const;
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc) const;
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const;
template <class FluidState>
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
template <class FluidState>
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
template <class FluidState>
void evalKrDeriv(const FluidState& /* fluidState */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
template <class FluidState>
void evalKrDeriv(const FluidState& /* fluidState */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
template <class FluidState>
void evalPc(const FluidState& /* fluidState */, double* /* pc */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
template <class FluidState>
void evalPcDeriv(const FluidState& /* fluidState */, double* /* pc */, double* /* dpcds */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
private:
};
typedef SatFuncStone2<UniformTableLinear<double> > SatFuncStone2Uniform;
typedef SatFuncStone2<NonuniformTableLinear<double> > SatFuncStone2Nonuniform;
template<class TableType>
template <class FluidState>
void SatFuncStone2<TableType>::evalKr(const FluidState& fluidState, double* kr) const
{
if (this->phase_usage.num_phases == 3) {
// Stone-II relative permeability model.
double sw = fluidState.saturation(BlackoilPhases::Aqua);
double sg = fluidState.saturation(BlackoilPhases::Vapour);
double krw = this->krw_(sw);
double krg = this->krg_(sg);
double krow = this->krow_(sw + sg); // = 1 - so
double krog = this->krog_(sg); // = 1 - so - sw
double krocw = this->krocw_;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[BlackoilPhases::Liquid] < 0.0) {
kr[BlackoilPhases::Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
template <class FluidState>
void SatFuncStone2<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Stone-II relative permeability model.
double sw = fluidState.saturation(BlackoilPhases::Aqua);
double sg = fluidState.saturation(BlackoilPhases::Vapour);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krow = this->krow_(sw + sg);
double dkrow = this->krow_.derivative(sw + sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
double krocw = this->krocw_;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[BlackoilPhases::Liquid] < 0.0) {
kr[BlackoilPhases::Liquid] = 0.0;
}
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
double dkrow = this->krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
template <class FluidState>
void SatFuncStone2<TableType>::evalPc(const FluidState& fluidState, double* pc) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(fluidState.saturation(pos));
}
}
template<class TableType>
template <class FluidState>
void SatFuncStone2<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = this->phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcow_.derivative(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcog_.derivative(fluidState.saturation(pos));
}
}
} // namespace Opm
#endif // SATFUNCSTONE2_HPP
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
-1
View File
@@ -23,7 +23,6 @@
#include <opm/core/props/satfunc/SaturationPropsInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SatFuncMultiplexer.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>