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

Merged upstream/master

Browse Source
This commit is contained in:
Joakim Hove 2014-02-04 21:42:07 +01:00
commit 6c8cba9d19
20 changed files with 2008 additions and 1824 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
cmake/Modules/opm-autodiff-prereqs.cmake
View File

@ -19,5 +19,5 @@ set (opm-autodiff_DEPS
dune-istl REQUIRED;
opm-core REQUIRED"
# Eigen
"Eigen3 3.1 REQUIRED"
"Eigen3 3.1.2 "
)

4
opm/core/grid/GridManager.cpp
View File

@ -216,6 +216,10 @@ namespace Opm
}
top_depths_vec.resize(x.size()*y.size(), tops[0]);
top_depths = &top_depths_vec[0];
} else {
OPM_THROW(std::runtime_error,
"Could not find either TOPS or DEPTHZ keyword, "
"one of them is required for initialization with DXV/DYV/DZV.");
}
// Construct grid.

29
opm/core/io/eclipse/EclipseGridParser.cpp
View File

@ -100,10 +100,17 @@ namespace EclipseKeywords
string("RV"),
string("DXV"), string("DYV"), string("DZV"),
string("DEPTHZ"), string("TOPS"), string("MAPAXES"),
string("SWCR"), string("SWL"), string("SWU"),
string("SOWCR"), string("KRW"), string("KRWR"),
string("KRO"), string("KRORW"), string("NTG"),
string("RHO")
string("SWL"), string("SWCR"), string("SWU"),
string("SGL"), string("SGCR"), string("SGU"),
string("SOWCR"), string("SOGCR"), string("KRW"),
string("KRWR"), string("KRG"), string("KRGR"),
string("KRO"), string("KRORW"), string("KRORG"),
string("ISWL"), string("ISWCR"), string("ISWU"),
string("ISGL"), string("ISGCR"), string("ISGU"),
string("ISOWCR"), string("ISOGCR"), string("IKRW"),
string("IKRWR"), string("IKRG"), string("IKRGR"),
string("IKRO"), string("IKRORW"), string("IKRORG"),
string("NTG"), string("RHO")
};
const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]);
@ -561,9 +568,17 @@ void EclipseGridParser::convertToSI()
key == "LAMEMOD" || key == "SHEARMOD" || key == "POISSONMOD" ||
key == "PWAVEMOD" || key == "MULTPV" || key == "PWAVEMOD" ||
key == "SGAS" || key == "SWAT" || key == "SOIL" ||
key == "NTG" || key == "SWCR" || key == "SWL" ||
key == "SWU" || key == "SOWCR" || key == "KRW" ||
key == "KRWR" || key == "KRORW" || key == "KRO" ||
key == "SWL" || key == "SWCR" || key == "SWU" ||
key == "SGL" || key == "SGCR" || key == "SGU" ||
key == "SOWCR" || key == "SOGCR" || key == "KRW" ||
key == "KRWR" || key == "KRG" || key == "KRGR" ||
key == "KRO" || key == "KRORW" || key == "KRORG" ||
key == "ISWL" || key == "ISWCR" || key == "ISWU" ||
key == "ISGL" || key == "ISGCR" || key == "ISGU" ||
key == "ISOWCR" || key == "ISOGCR" || key == "IKRW" ||
key == "IKRWR" || key == "IKRG" || key == "IKRGR" ||
key == "IKRO" || key == "IKRORW" || key == "IKRORG" ||
key == "NTG" ||
key == "RHO") /* nonstandard field with no unit logic. use with caution */ {
unit = 1.0;
do_convert = false; // Dimensionless keywords...

92
opm/core/io/eclipse/EclipseWriter.cpp
View File

@ -575,6 +575,12 @@ private:
fortio_fclose) { }
};
} // anonymous namespace
// Note: the following parts were taken out of the anonymous
// namespace, since EclipseSummary is now used as a pointer member in
// EclipseWriter and forward declared in EclipseWriter.hpp.
// forward decl. of mutually dependent type
struct EclipseWellReport;
@ -655,6 +661,7 @@ private:
}
};
// in order to get RTTI for this "class" (which is just a typedef), we must
// ask the compiler to explicitly instantiate it.
template struct EclipseHandle<ecl_sum_tstep_struct>;
@ -670,7 +677,7 @@ protected:
PhaseUsage uses, /* phases present */
BlackoilPhases::PhaseIndex phase, /* oil, water or gas */
WellType type, /* prod. or inj. */
char aggregation, /* rate or total */
char aggregation, /* rate or total or BHP */
std::string unit)
: EclipseHandle <smspec_node_type> (
ecl_sum_add_var (summary,
@ -718,22 +725,26 @@ private:
char aggregation) {
std::string name;
name += 'W'; // well
switch (phase) {
if (aggregation == 'B') {
name += "BHP";
} else {
switch (phase) {
case BlackoilPhases::Aqua: name += 'W'; break; /* water */
case BlackoilPhases::Vapour: name += 'G'; break; /* gas */
case BlackoilPhases::Liquid: name += 'O'; break; /* oil */
default:
OPM_THROW(std::runtime_error,
"Unknown phase used in blackoil reporting");
}
switch (type) {
}
switch (type) {
case WellType::INJECTOR: name += 'I'; break;
case WellType::PRODUCER: name += 'P'; break;
default:
OPM_THROW(std::runtime_error,
"Unknown well type used in blackoil reporting");
}
name += aggregation; /* rate ('R') or total ('T') */
}
name += aggregation; /* rate ('R') or total ('T') */
return name;
}
protected:
@ -743,6 +754,13 @@ protected:
const double value = sign_ * wellState.wellRates () [index_] * convFactor;
return value;
}
double bhp (const WellState& wellstate) {
// Note that 'index_' is used here even though it is meant
// to give a (well,phase) pair.
const int num_phases = wellstate.wellRates().size() / wellstate.bhp().size();
return wellstate.bhp()[index_/num_phases];
}
};
/// Monitors the rate given by a well.
@ -805,6 +823,31 @@ private:
double total_;
};
/// Monitors the bottom hole pressure in a well.
struct EclipseWellBhp : public EclipseWellReport {
EclipseWellBhp (const EclipseSummary& summary,
const EclipseGridParser& parser,
int whichWell,
PhaseUsage uses,
BlackoilPhases::PhaseIndex phase,
WellType type)
: EclipseWellReport (summary,
parser,
whichWell,
uses,
phase,
type,
'B',
"Pascal")
{ }
virtual double update (const SimulatorTimer& timer,
const WellState& wellState)
{
return bhp(wellState);
}
};
inline void
EclipseSummary::writeTimeStep (const SimulatorTimer& timer,
const WellState& wellState) {
@ -861,8 +904,30 @@ EclipseSummary::addWells (const EclipseGridParser& parser,
}
}
}
// Add BHP monitors
for (int whichWell = 0; whichWell != numWells; ++whichWell) {
// In the call below: uses, phase and the well type arguments
// are not used, except to set up an index that stores the
// well indirectly. For details see the implementation of the
// EclipseWellReport constructor, and the method
// EclipseWellReport::bhp().
BlackoilPhases::PhaseIndex phase = BlackoilPhases::Liquid;
if (!uses.phase_used[BlackoilPhases::Liquid]) {
phase = BlackoilPhases::Vapour;
}
add (std::unique_ptr <EclipseWellReport> (
new EclipseWellBhp (*this,
parser,
whichWell,
uses,
phase,
WELL_TYPES[0])));
}
}
namespace {
/// Helper method that can be used in keyword transformation (must curry
/// the barsa argument)
static double toBar (const double& pressure) {
@ -902,6 +967,11 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer,
/* Initial solution (pressure and saturation) */
writeSolution (timer, reservoirState, wellState);
/* Create summary object (could not do it at construction time,
since it requires knowledge of the start time). */
summary_.reset(new EclipseSummary(outputDir_, baseName_, timer, *parser_));
summary_->addWells (*parser_, uses_);
}
void EclipseWriter::writeSolution (const SimulatorTimer& timer,
@ -952,9 +1022,15 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
// (first timestep, in practice), and reused later. but how to do this
// without keeping the complete summary in memory (which will then
// accumulate all the timesteps)?
EclipseSummary sum (outputDir_, baseName_, timer, *parser_);
sum.addWells (*parser_, uses_);
sum.writeTimeStep (timer, wellState);
//
// Note: The answer to the question above is still not settled,
// but now we do keep the complete summary in memory, as a member
// variable in the EclipseWriter class, instead of creating a
// temporary EclipseSummary in this function every time it is
// called. This has been changed so that the final summary file
// will contain data from the whole simulation, instead of just
// the last step.
summary_->writeTimeStep (timer, wellState);
}
#else

2
opm/core/io/eclipse/EclipseWriter.hpp
View File

@ -28,6 +28,7 @@
#include <memory> // std::unique_ptr
struct UnstructuredGrid;
struct EclipseSummary;
namespace Opm {
@ -90,6 +91,7 @@ private:
std::string outputDir_;
std::string baseName_;
PhaseUsage uses_; // active phases in the input deck
std::shared_ptr <EclipseSummary> summary_;
/// Write solution field variables (pressure and saturation)
void writeSolution (const SimulatorTimer& timer,

39
opm/core/io/eclipse/SpecialEclipseFields.hpp
View File

@ -2313,8 +2313,8 @@ struct ENPTVD : public SpecialBase {
}
virtual void read(std::istream & is) {
table_.resize(5);
std::vector<std::vector<double> > sub_table(5);
table_.resize(9);
std::vector<std::vector<double> > sub_table(9);
while (!is.eof()) {
if(is.peek() == int('/')) {
if (sub_table[0].empty() && !(table_[0].empty())) {
@ -2332,19 +2332,14 @@ struct ENPTVD : public SpecialBase {
if (data[0] == -1.0) {
OPM_THROW(std::runtime_error, "Error reading ENPTVD data - depth can not be defaulted.");
}
if ((data[4] != -1.0) || (data[5] != -1.0) || (data[6] != -1.0) || (data[8] != -1.0)) {
OPM_THROW(std::runtime_error, "Error reading ENPTVD data - non-default values in column 5-7,9 not supported.");
for(std::vector<std::vector<double> >::size_type i=0; i<sub_table.size(); ++i) {
sub_table[i].push_back(data[i]); // [0-8]: depth swl swcr swu sgl sgcr sgu sowcr sogcr
}
sub_table[0].push_back(data[0]); //depth
sub_table[1].push_back(data[1]); //swl
sub_table[2].push_back(data[2]); //swcr
sub_table[3].push_back(data[3]); //swu
sub_table[4].push_back(data[7]); //sowcr
is >> ignoreWhitespace;
if(is.peek() == int('/')) {
is >> ignoreLine;
if (sub_table[0].size() >= 2) {
insertDefaultValues(sub_table, 5, -1.0, false);
insertDefaultValues(sub_table, 9, -1.0, false);
std::vector<std::vector<double> >::iterator it_sub = sub_table.begin();
for(std::vector<std::vector<std::vector<double> > >::size_type i=0; i<table_.size(); ++i) {
table_[i].push_back(*it_sub);
@ -2369,9 +2364,9 @@ struct ENPTVD : public SpecialBase {
virtual void write(std::ostream & os) const {
os << name() << '\n';
std::cout << "-----depth-------swl------swcr-------swu-----sowcr" << std::endl;
std::cout << "-----depth-------swl------swcr-------swu-------sgl------sgcr-------sgu-----sowcr-----sogcr" << std::endl;
for (std::vector<std::vector<double> >::size_type j=0; j<table_[0].size(); ++j) {
std::cout << "--------------------------------------------------" << std::endl;
std::cout << "------------------------------------------------------------------------------------------" << std::endl;
for (std::vector<double>::size_type k=0; k<table_[0][j].size(); ++k) {
for (std::vector<std::vector<std::vector<double> > >::size_type i=0; i<table_.size(); ++i) {
std::cout << std::setw(10) << table_[i][j][k];
@ -2400,8 +2395,8 @@ struct ENKRVD : public SpecialBase {
}
virtual void read(std::istream & is) {
table_.resize(5);
std::vector<std::vector<double> > sub_table(5);
table_.resize(8);
std::vector<std::vector<double> > sub_table(8);
while (!is.eof()) {
if(is.peek() == int('/')) {
if (sub_table[0].empty() && !(table_[0].empty())) {
@ -2419,19 +2414,15 @@ struct ENKRVD : public SpecialBase {
if (data[0] == -1.0) {
OPM_THROW(std::runtime_error, "Error reading ENKRVD data - depth can not be defaulted.");
}
if ((data[2] != -1.0) || (data[5] != -1.0) || (data[6] != -1.0)) {
OPM_THROW(std::runtime_error, "Error reading ENKRVD data - non-default values in column 3,6-7 not supported.");
for(std::vector<std::vector<double> >::size_type i=0; i<sub_table.size(); ++i) {
sub_table[i].push_back(data[i]); // [0-7]: depth krw krg kro krwr krgr krorw krorg
}
sub_table[0].push_back(data[0]); //depth
sub_table[1].push_back(data[1]); //krw
sub_table[2].push_back(data[3]); //kro
sub_table[3].push_back(data[4]); //krw(sowcr)
sub_table[4].push_back(data[7]); //kro(swcr)
is >> ignoreWhitespace;
if(is.peek() == int('/')) {
is >> ignoreLine;
if (sub_table[0].size() >= 2) {
insertDefaultValues(sub_table, 5, -1.0, false);
insertDefaultValues(sub_table, 8, -1.0, false);
std::vector<std::vector<double> >::iterator it_sub = sub_table.begin();
for(std::vector<std::vector<std::vector<double> > >::size_type i=0; i<table_.size(); ++i) {
table_[i].push_back(*it_sub);
@ -2457,9 +2448,9 @@ struct ENKRVD : public SpecialBase {
virtual void write(std::ostream & os) const {
os << name() << '\n';
std::cout << "-----depth-------krw------krow------krwr-----krorw" << std::endl;
std::cout << "-----depth-------krw-------krg-------kro------krwr------krgr-----krorw-----krorg" << std::endl;
for (std::vector<std::vector<double> >::size_type j=0; j<table_[0].size(); ++j) {
std::cout << "--------------------------------------------------" << std::endl;
std::cout << "--------------------------------------------------------------------------------" << std::endl;
for (std::vector<double>::size_type k=0; k<table_[0][j].size(); ++k) {
for (std::vector<std::vector<std::vector<double> > >::size_type i=0; i<table_.size(); ++i) {
std::cout << std::setw(10) << table_[i][j][k];

312
opm/core/props/satfunc/SatFuncBase.hpp Normal file
View File

@ -0,0 +1,312 @@
/*
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/io/eclipse/EclipseGridParser.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 <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))
void printMe(std::ostream & out);
};
Transform wat;
Transform watoil;
Transform gas;
Transform gasoil;
};
// Hysteresis
struct SatHyst {
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:
void init(const EclipseGridParser& deck,
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.
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(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
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;
}
}
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// 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;
}
}
}
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

503
opm/core/props/satfunc/SatFuncGwseg.cpp
View File

@ -29,504 +29,5 @@
namespace Opm
{
void SatFuncGwsegUniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
// 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);
SatFuncGwsegUniform::ExtendTable(sw,sw_ex,1);
SatFuncGwsegUniform::ExtendTable(krw,krw_ex,0);
SatFuncGwsegUniform::ExtendTable(krow,krow_ex,0);
SatFuncGwsegUniform::ExtendTable(pcow,pcow_ex,0);
buildUniformMonotoneTable(sw_ex, krw_ex, samples, krw_);
buildUniformMonotoneTable(sw_ex, krow_ex, samples, krow_);
buildUniformMonotoneTable(sw_ex, pcow_ex, samples, pcow_);
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();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// 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);
SatFuncGwsegUniform::ExtendTable(sg,sg_ex,1);
SatFuncGwsegUniform::ExtendTable(krg,krg_ex,0);
SatFuncGwsegUniform::ExtendTable(krog,krog_ex,0);
SatFuncGwsegUniform::ExtendTable(pcog,pcog_ex,0);
buildUniformMonotoneTable(sg_ex, krg_ex, samples, krg_);
buildUniformMonotoneTable(sg_ex, krog_ex, samples, krog_);
buildUniformMonotoneTable(sg_ex, pcog_ex, samples, pcog_);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncGwsegUniform::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];
}
}
void SatFuncGwsegUniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
const double swco = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[Vapour];
// xw and xg are the fractions occupied by water and gas zones.
const double eps = 1e-6;
const double xw = (sw - swco) / std::max(sg + sw - swco, eps);
const double xg = 1 - xw;
const double ssw = sg + sw;
const double ssg = sw - swco + sg;
const double krw = krw_(ssw);
const double krg = krg_(ssg);
const double krow = krow_(ssw);
const double krog = krog_(ssg);
kr[Aqua] = xw*krw;
kr[Vapour] = xg*krg;
kr[Liquid] = xw*krow + xg*krog;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncGwsegUniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = 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.
const double swco = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[Vapour];
// xw and xg are the fractions occupied by water and gas zones.
const double eps = 1e-6;
const double ssw = sg + sw;
const double ssg = std::max(sg + sw - swco, eps);
const double krw = krw_(ssw);
const double krg = krg_(ssg);
const double krow = krow_(ssw);
const double krog = krog_(ssg);
const double xw = (sw - swco) / ssg;
const double xg = 1 - xw;
kr[Aqua] = xw*krw;
kr[Vapour] = xg*krg;
kr[Liquid] = xw*krow + xg*krog;
// Derivatives.
const double dkrww = krw_.derivative(ssw);
const double dkrgg = krg_.derivative(ssg);
const double dkrow = krow_.derivative(ssw);
const double dkrog = krog_.derivative(ssg);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
dkrds[Aqua + Aqua*np] = (xg/d)*krw + xw*dkrww;
dkrds[Aqua + Vapour*np] = -(xw/d)*krw + xw*dkrww;
dkrds[Liquid + Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[Liquid + Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[Vapour + Aqua*np] = -(xg/d)*krg + xg*dkrgg;
dkrds[Vapour + Vapour*np] = (xw/d)*krg + xg*dkrgg;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = 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(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncGwsegUniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncGwsegUniform::evalPcDeriv(const double* s, 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 = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
// ====== Methods for SatFuncGwsegNonuniform ======
void SatFuncGwsegNonuniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int /*samples*/)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
// 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);
SatFuncGwsegNonuniform::ExtendTable(sw,sw_ex,1);
SatFuncGwsegNonuniform::ExtendTable(krw,krw_ex,0);
SatFuncGwsegNonuniform::ExtendTable(krow,krow_ex,0);
SatFuncGwsegNonuniform::ExtendTable(pcow,pcow_ex,0);
krw_ = NonuniformTableLinear<double>(sw_ex, krw_ex);
krow_ = NonuniformTableLinear<double>(sw_ex, krow_ex);
pcow_ = NonuniformTableLinear<double>(sw_ex, pcow_ex);
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();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// 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);
SatFuncGwsegNonuniform::ExtendTable(sg,sg_ex,1);
SatFuncGwsegNonuniform::ExtendTable(krg,krg_ex,0);
SatFuncGwsegNonuniform::ExtendTable(krog,krog_ex,0);
SatFuncGwsegNonuniform::ExtendTable(pcog,pcog_ex,0);
krg_ = NonuniformTableLinear<double>(sg_ex, krg_ex);
krog_ = NonuniformTableLinear<double>(sg_ex, krog_ex);
pcog_ = NonuniformTableLinear<double>(sg_ex, pcog_ex);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncGwsegNonuniform::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];
}
}
void SatFuncGwsegNonuniform::evalKr(const double* s, double* kr) const
{
if (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 = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[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 = krw_(sw);
const double krg = krg_(sg);
const double krow = krow_(ssw);
const double krog = krog_(ssg);
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = xw*krow + xg*krog;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncGwsegNonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = 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 = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[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 = krw_(sw);
const double krg = krg_(sg);
const double krow = krow_(ssw);
const double krog = krog_(d);
const double xg = 1 - xw;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = xw*krow + xg*krog;
// Derivatives.
const double dkrww = krw_.derivative(sw);
const double dkrgg = krg_.derivative(sg);
const double dkrow = krow_.derivative(ssw);
const double dkrog = krog_.derivative(d);
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Liquid + Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[Liquid + Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[Vapour + Vapour*np] = dkrgg;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = 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(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncGwsegNonuniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncGwsegNonuniform::evalPcDeriv(const double* s, 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 = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
} // namespace Opm
// SatFuncGwseg.cpp can now be removed
} // namespace Opm

643
opm/core/props/satfunc/SatFuncGwseg.hpp
View File

@ -19,81 +19,602 @@
#ifndef SATFUNCGWSEG_HPP
#define SATFUNCGWSEG_HPP
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <vector>
#include <opm/core/props/satfunc/SatFuncBase.hpp>
namespace Opm
{
class SatFuncGwsegUniform : public BlackoilPhases
template<class TableType>
class SatFuncGwseg : public SatFuncBase<TableType>
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
void ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
void evalKr(const double* s, double* kr, const EPSTransforms* epst) const;
void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) 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:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
UniformTableLinear<double> krw_;
UniformTableLinear<double> krow_;
UniformTableLinear<double> pcow_;
UniformTableLinear<double> krg_;
UniformTableLinear<double> krog_;
UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
typedef SatFuncGwseg<UniformTableLinear<double> > SatFuncGwsegUniform;
typedef SatFuncGwseg<NonuniformTableLinear<double> > SatFuncGwsegNonuniform;
class SatFuncGwsegNonuniform : public BlackoilPhases
template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, double* kr) const
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
void ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) 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(s[BlackoilPhases::Aqua], swco);
const double sg = s[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 = s[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 = s[gpos];
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, 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(s[BlackoilPhases::Aqua], swco);
const double sg = s[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 = s[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 = s[gpos];
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, 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(s[BlackoilPhases::Aqua], swco);
const double sg = s[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 = s[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 = s[gpos];
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, 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(s[BlackoilPhases::Aqua], swco);
const double sg = s[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 = s[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 = s[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>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, 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(s[BlackoilPhases::Aqua], swco);
const double sg = s[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_(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;
// 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(1.0-_ssow));
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(1.0-_ssog-_swco));
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 = s[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 = s[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>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, 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(s[BlackoilPhases::Aqua], swco);
const double sg = s[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_(1.0-_ssow), this->krorw_);
dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(1.0-_ssow));
} 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_(1.0-_ssow), this->krorw_);
dkrow = _dsdssow*epst_hyst->watoil.scaleKrDeriv(ssow_shifted, this->krow_.derivative(1.0-_ssow));
}
// 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_(1.0-_ssog-_swco), this->krorg_);
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(1.0-_ssog-_swco));
// 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 = s[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 = s[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>
void SatFuncGwseg<TableType>::evalPc(const double* s, 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_(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPc(const double* s, 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(s[pos], this->smin_[pos], this->smax_[pos]);
pc[pos] = 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(s[pos], this->smin_[pos], this->smax_[pos]);
pc[pos] = this->pcog_(_sg);
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPcDeriv(const double* s, 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_(s[pos]);
dpcds[np*pos + pos] = this->pcow_.derivative(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]);
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPcDeriv(const double* s, 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(s[pos], this->smin_[pos], this->smax_[pos]);
pc[pos] = this->pcow_(s[pos]);
double _dsdsw = epst->wat.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]);
dpcds[np*pos + pos] = _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(s[pos], this->smin_[pos], this->smax_[pos]);
pc[pos] = this->pcog_(_sg);
double _dsdsg = epst->gas.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]);
dpcds[np*pos + pos] = _dsdsg*this->pcog_.derivative(_sg);
}
}
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
NonuniformTableLinear<double> krw_;
NonuniformTableLinear<double> krow_;
NonuniformTableLinear<double> pcow_;
NonuniformTableLinear<double> krg_;
NonuniformTableLinear<double> krog_;
NonuniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
} // namespace Opm
#endif // SATFUNCGWSEG_HPP

646
opm/core/props/satfunc/SatFuncSimple.cpp
View File

@ -30,536 +30,194 @@
namespace Opm
{
void SatFuncSimpleUniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
// 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)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
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);
SatFuncSimpleUniform::ExtendTable(sw,sw_ex,1);
SatFuncSimpleUniform::ExtendTable(krw,krw_ex,0);
SatFuncSimpleUniform::ExtendTable(krow,krow_ex,0);
SatFuncSimpleUniform::ExtendTable(pcow,pcow_ex,0);
buildUniformMonotoneTable(sw_ex, krw_ex, samples, krw_);
buildUniformMonotoneTable(sw_ex, krow_ex, samples, krow_);
buildUniformMonotoneTable(sw_ex, pcow_ex, samples, pcow_);
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;
}
}
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// 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);
SatFuncSimpleUniform::ExtendTable(sg,sg_ex,1);
SatFuncSimpleUniform::ExtendTable(krg,krg_ex,0);
SatFuncSimpleUniform::ExtendTable(krog,krog_ex,0);
SatFuncSimpleUniform::ExtendTable(pcog,pcog_ex,0);
buildUniformMonotoneTable(sg_ex, krg_ex, samples, krg_);
buildUniformMonotoneTable(sg_ex, krog_ex, samples, krog_);
buildUniformMonotoneTable(sg_ex, pcog_ex, samples, pcog_);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
table = NonuniformTableLinear<double>(arg, value);
}
void SatFuncSimpleUniform::ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const
template<>
void SatFuncBase<UniformTableLinear<double> >::initializeTableType(UniformTableLinear<double> & table,
const std::vector<double>& arg,
const std::vector<double>& value,
const int samples)
{
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];
}
buildUniformMonotoneTable(arg, value, samples, table);
}
void SatFuncSimpleUniform::evalKr(const double* s, double* kr) const
double EPSTransforms::Transform::scaleSat(double s, double s_r, double s_cr, double s_max) const
{
if (phase_usage.num_phases == 3) {
// A simplified relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double krg = krg_(sg);
double krow = krow_(sw + sg); // = 1 - so
// double krog = krog_(sg); // = 1 - so - sw
// double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krow;
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
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;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double so = s[opos];
double krow = krow_(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
} else if (s <= sr) {
return (s <= scr) ? s_cr : s_cr+(s-scr)*slope1;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
return (s >= smax) ? s_max : s_r+(s-sr)*slope2;
}
}
void SatFuncSimpleUniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
double EPSTransforms::Transform::scaleSatInv(double ss, double s_r, double s_cr, double s_max) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// A simplified relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krow = krow_(sw + sg);
double dkrow = krow_.derivative(sw + sg);
// double krog = krog_(sg);
// double dkrog = krog_.derivative(sg);
// double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krow;
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
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;
}
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Vapour + Vapour*np] = dkrgg;
//dkrds[Liquid + Aqua*np] = dkrow;
dkrds[Liquid + Liquid*np] = -dkrow;
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[Liquid + 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 (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double so = s[opos];
double krow = krow_(1.0-so);
double dkrow = 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 if (ss <= s_r) {
return (ss <= s_cr) ? scr : scr+(ss-s_cr)/slope1;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncSimpleUniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
return (ss >= s_max) ? smax : sr+(ss-s_r)/slope2;
}
}
void SatFuncSimpleUniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
double EPSTransforms::Transform::scaleSatDeriv(double s, double s_r, double s_cr, double s_max) 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 = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
// ====== Methods for SatFuncSimpleNonuniform ======
void SatFuncSimpleNonuniform::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];
}
}
void SatFuncSimpleNonuniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int /*samples*/)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
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)");
if (doNotScale) {
return 1.0;
} else if (!do_3pt) { // 2-pt
if (s <= scr) {
return 0.0;
} else {
return (s >= smax) ? 0.0 : slope1;
}
// 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);
SatFuncSimpleNonuniform::ExtendTable(sw,sw_ex,1);
SatFuncSimpleNonuniform::ExtendTable(krw,krw_ex,0);
SatFuncSimpleNonuniform::ExtendTable(krow,krow_ex,0);
SatFuncSimpleNonuniform::ExtendTable(pcow,pcow_ex,0);
krw_ = NonuniformTableLinear<double>(sw_ex, krw_ex);
krow_ = NonuniformTableLinear<double>(sw_ex, krow_ex);
pcow_ = NonuniformTableLinear<double>(sw_ex, pcow_ex);
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;
}
}
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// 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);
SatFuncSimpleNonuniform::ExtendTable(sg,sg_ex,1);
SatFuncSimpleNonuniform::ExtendTable(krg,krg_ex,0);
SatFuncSimpleNonuniform::ExtendTable(krog,krog_ex,0);
SatFuncSimpleNonuniform::ExtendTable(pcog,pcog_ex,0);
krg_ = NonuniformTableLinear<double>(sg_ex, krg_ex);
krog_ = NonuniformTableLinear<double>(sg_ex, krog_ex);
pcog_ = NonuniformTableLinear<double>(sg_ex, pcog_ex);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncSimpleNonuniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// A simplified relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double krg = krg_(sg);
double krow = krow_(sw + sg); // = 1 - so
// double krog = krog_(sg); // = 1 - so - sw
// double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krow;
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double so = s[opos];
double krow = krow_(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
} else if (s <= sr) {
return (s <= scr) ? 0.0 : slope1;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
return (s >= smax) ? 0.0 : slope2;
}
}
void SatFuncSimpleNonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
double EPSTransforms::Transform::scaleSatPc(double s, double s_min, double s_max) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// A simplified relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krow = krow_(sw + sg);
double dkrow = krow_.derivative(sw + sg);
// double krog = krog_(sg);
// double dkrog = krog_.derivative(sg);
// double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krow;
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Vapour + Vapour*np] = dkrgg;
//dkrds[Liquid + Aqua*np] = dkrow;
dkrds[Liquid + Liquid*np] = -dkrow;
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[Liquid + 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 (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double so = s[opos];
double krow = krow_(1.0-so);
double dkrow = 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.
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 {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
return s_max;
}
}
void SatFuncSimpleNonuniform::evalPc(const double* s, double* pc) const
double EPSTransforms::Transform::scaleSatDerivPc(double s, double s_min, double s_max) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
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;
}
}
void SatFuncSimpleNonuniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
double EPSTransforms::Transform::scaleKr(double s, double kr, double krsr_tab) 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 = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
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

306
opm/core/props/satfunc/SatFuncSimple.hpp
View File

@ -19,80 +19,264 @@
#ifndef SATFUNCSIMPLE_HPP
#define SATFUNCSIMPLE_HPP
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <vector>
#include <opm/core/props/satfunc/SatFuncBase.hpp>
namespace Opm
{
class SatFuncSimpleUniform : public BlackoilPhases
template<class TableType>
class SatFuncSimple : public SatFuncBase<TableType>
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
void ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
void evalKr(const double* s, double* kr, const EPSTransforms* epst) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const;
void evalKrDeriv(const double* s, 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 ...");}
void evalPc(const double* s, double* pc, const EPSTransforms* epst) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
UniformTableLinear<double> krw_;
UniformTableLinear<double> krow_;
UniformTableLinear<double> pcow_;
UniformTableLinear<double> krg_;
UniformTableLinear<double> krog_;
UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
typedef SatFuncSimple<UniformTableLinear<double> > SatFuncSimpleUniform;
typedef SatFuncSimple<NonuniformTableLinear<double> > SatFuncSimpleNonuniform;
class SatFuncSimpleNonuniform : public BlackoilPhases
template<class TableType>
void SatFuncSimple<TableType>::evalKr(const double* s, double* kr) const
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
void ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
NonuniformTableLinear<double> krw_;
NonuniformTableLinear<double> krow_;
NonuniformTableLinear<double> pcow_;
NonuniformTableLinear<double> krg_;
NonuniformTableLinear<double> krog_;
NonuniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
if (this->phase_usage.num_phases == 3) {
// A simplified relative permeability model.
double sw = s[BlackoilPhases::Aqua];
double sg = s[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 = s[wpos];
double krw = this->krw_(sw);
double so = s[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 = s[gpos];
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
void SatFuncSimple<TableType>::evalKrDeriv(const double* s, 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 = s[BlackoilPhases::Aqua];
double sg = s[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 = s[wpos];
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double so = s[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 = s[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>
void SatFuncSimple<TableType>::evalKrDeriv(const double* s, 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(s[wpos], 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _dsdsw = epst->wat.scaleSatDeriv(s[wpos], 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _sg = epst->gas.scaleSat(s[gpos], 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _dsdsg = epst->gas.scaleSatDeriv(s[gpos], 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _sow = epst->watoil.scaleSat(1.0-s[wpos]-s[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-s[wpos]-s[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-s[wpos]-s[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-s[wpos]-s[gpos], 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double krw = epst->wat.scaleKr(s[wpos], this->krw_(_sw), this->krwr_);
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(s[wpos], this->krw_.derivative(_sw));
double krg = epst->gas.scaleKr(s[gpos], this->krg_(_sg), this->krgr_);
double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(s[gpos], this->krg_.derivative(_sg));
// TODO Check the arguments to the krow- and krog-tables below...
double krow = epst->watoil.scaleKr(1.0-s[wpos]-s[gpos], this->krow_(1.0-_sow-this->smin_[gpos]), this->krorw_); // ????
double dkrow = _dsdsow*epst->watoil.scaleKrDeriv(1.0-s[wpos]-s[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-s[wpos]-s[gpos], this->krorg_); // ????
//double dkrog = _dsdsog*epst->gasoil.scaleKrDeriv(1.0-s[wpos]-s[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 = s[wpos];
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double so = s[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 = s[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>
void SatFuncSimple<TableType>::evalPc(const double* s, 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_(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
}
}
template<class TableType>
void SatFuncSimple<TableType>::evalPcDeriv(const double* s, 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_(s[pos]);
dpcds[np*pos + pos] = this->pcow_.derivative(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]);
}
}
} // namespace Opm
#endif // SATFUNCSIMPLE_HPP

389
opm/core/props/satfunc/SatFuncStone2.cpp
View File

@ -29,390 +29,5 @@
namespace Opm
{
void SatFuncStone2Uniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
buildUniformMonotoneTable(sw, krw, samples, krw_);
buildUniformMonotoneTable(sw, krow, samples, krow_);
buildUniformMonotoneTable(sw, pcow, samples, pcow_);
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();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
buildUniformMonotoneTable(sg, krg, samples, krg_);
buildUniformMonotoneTable(sg, krog, samples, krog_);
buildUniformMonotoneTable(sg, pcog, samples, pcog_);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncStone2Uniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double krg = krg_(sg);
double krow = krow_(sw + sg); // = 1 - so
double krog = krog_(sg); // = 1 - so - sw
double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncStone2Uniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krow = krow_(sw + sg);
double dkrow = krow_.derivative(sw + sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Vapour + Vapour*np] = dkrgg;
dkrds[Liquid + Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[Liquid + 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 (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = 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(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncStone2Uniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncStone2Uniform::evalPcDeriv(const double* s, 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 = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
// ====== Methods for SatFuncStone2Nonuniform ======
void SatFuncStone2Nonuniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int /*samples*/)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
krw_ = NonuniformTableLinear<double>(sw, krw);
krow_ = NonuniformTableLinear<double>(sw, krow);
pcow_ = NonuniformTableLinear<double>(sw, pcow);
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();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
krg_ = NonuniformTableLinear<double>(sg, krg);
krog_ = NonuniformTableLinear<double>(sg, krog);
pcog_ = NonuniformTableLinear<double>(sg, pcog);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncStone2Nonuniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double krg = krg_(sg);
double krow = krow_(sw + sg); // = 1 - so
double krog = krog_(sg); // = 1 - so - sw
double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncStone2Nonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krow = krow_(sw + sg);
double dkrow = krow_.derivative(sw + sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Vapour + Vapour*np] = dkrgg;
dkrds[Liquid + Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[Liquid + 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 (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = 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(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncStone2Nonuniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncStone2Nonuniform::evalPcDeriv(const double* s, 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 = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
} // namespace Opm
// SatFuncStone2.cpp can now be removed
} // namespace Opm

216
opm/core/props/satfunc/SatFuncStone2.hpp
View File

@ -19,74 +19,180 @@
#ifndef SATFUNCSTONE2_HPP
#define SATFUNCSTONE2_HPP
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <vector>
#include <opm/core/props/satfunc/SatFuncBase.hpp>
namespace Opm
{
class SatFuncStone2Uniform : public BlackoilPhases
template<class TableType>
class SatFuncStone2 : public SatFuncBase<TableType>
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
void evalKr(const double* s, double* kr, const EPSTransforms* epst) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalKrDeriv(const double* s, 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 ...");}
void evalPc(const double* s, double* pc, const EPSTransforms* epst) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
UniformTableLinear<double> krw_;
UniformTableLinear<double> krow_;
UniformTableLinear<double> pcow_;
UniformTableLinear<double> krg_;
UniformTableLinear<double> krog_;
UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
typedef SatFuncStone2<UniformTableLinear<double> > SatFuncStone2Uniform;
typedef SatFuncStone2<NonuniformTableLinear<double> > SatFuncStone2Nonuniform;
class SatFuncStone2Nonuniform : public BlackoilPhases
template<class TableType>
void SatFuncStone2<TableType>::evalKr(const double* s, double* kr) const
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
NonuniformTableLinear<double> krw_;
NonuniformTableLinear<double> krow_;
NonuniformTableLinear<double> pcow_;
NonuniformTableLinear<double> krg_;
NonuniformTableLinear<double> krog_;
NonuniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
if (this->phase_usage.num_phases == 3) {
// Stone-II relative permeability model.
double sw = s[BlackoilPhases::Aqua];
double sg = s[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 = s[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 = s[gpos];
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
void SatFuncStone2<TableType>::evalKrDeriv(const double* s, 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 = s[BlackoilPhases::Aqua];
double sg = s[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 = s[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 = s[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>
void SatFuncStone2<TableType>::evalPc(const double* s, 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_(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
}
}
template<class TableType>
void SatFuncStone2<TableType>::evalPcDeriv(const double* s, 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_(s[pos]);
dpcds[np*pos + pos] = this->pcow_.derivative(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]);
}
}
} // namespace Opm
#endif // SATFUNCSTONE2_HPP

50
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
View File

@ -101,34 +101,54 @@ namespace Opm
const int* cells,
double* smin,
double* smax) const;
/// Update saturation state for the hysteresis tracking
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
void updateSatHyst(const int n,
const int* cells,
const double* s);
private:
PhaseUsage phase_usage_;
std::vector<SatFuncSet> satfuncset_;
std::vector<int> cell_to_func_; // = SATNUM - 1
struct { // End point scaling parameters
std::vector<double> swl_;
std::vector<double> swcr_;
std::vector<double> swu_;
std::vector<double> sowcr_;
std::vector<double> krw_;
std::vector<double> krwr_;
std::vector<double> kro_;
std::vector<double> krorw_;
} eps_;
bool do_eps_; // ENDSCALE is active
bool do_3pt_; // SCALECRS: YES~true NO~false
bool do_hyst_; // Keywords ISWL etc detected
std::vector<EPSTransforms> eps_transf_;
std::vector<EPSTransforms> eps_transf_hyst_;
std::vector<SatHyst> sat_hyst_;
typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
void initEPS(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const std::string& keyword,
std::vector<double>& scaleparam);
void relpermEPS(const double *s, const int cell, double *kr, double *dkrds= 0) const;
const UnstructuredGrid& grid);
void initEPSHyst(const EclipseGridParser& deck,
const UnstructuredGrid& grid);
void initEPSKey(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const std::string& keyword,
std::vector<double>& scaleparam);
void initEPSParam(const int cell,
EPSTransforms::Transform& data,
const bool oil,
const double sl_tab,
const double scr_tab,
const double su_tab,
const double sxcr_tab,
const double s0_tab,
const double krsr_tab,
const double krmax_tab,
const std::vector<double>& sl,
const std::vector<double>& scr,
const std::vector<double>& su,
const std::vector<double>& sxcr,
const std::vector<double>& s0,
const std::vector<double>& krsr,
const std::vector<double>& krmax);
};

581
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
View File

@ -103,9 +103,9 @@ namespace Opm
do_eps_ = false;
do_3pt_ = false;
if (deck.hasField("ENDSCALE")) {
if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) {
OPM_THROW(std::runtime_error, "Currently endpoint-scaling limited to oil-water systems without gas.");
}
//if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) {
// OPM_THROW(std::runtime_error, "Currently endpoint-scaling limited to oil-water systems without gas.");
//}
if (deck.getENDSCALE().dir_switch_ != std::string("NODIR")) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'NODIR' accepted.");
}
@ -118,14 +118,23 @@ namespace Opm
}
}
do_eps_ = true;
initEPS(deck, grid, std::string("SWCR"), eps_.swcr_);
initEPS(deck, grid, std::string("SWL"), eps_.swl_);
initEPS(deck, grid, std::string("SWU"), eps_.swu_);
initEPS(deck, grid, std::string("SOWCR"), eps_.sowcr_);
initEPS(deck, grid, std::string("KRW"), eps_.krw_);
initEPS(deck, grid, std::string("KRWR"), eps_.krwr_);
initEPS(deck, grid, std::string("KRO"), eps_.kro_);
initEPS(deck, grid, std::string("KRORW"), eps_.krorw_);
initEPS(deck, grid);
// For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR
do_hyst_ = deck.hasField("ISWL") || deck.hasField("ISWU") || deck.hasField("ISWCR") || deck.hasField("ISGL") ||
deck.hasField("ISGU") || deck.hasField("ISGCR") || deck.hasField("ISOWCR") || deck.hasField("ISOGCR");
if (do_hyst_) {
if (deck.hasField("KRW") || deck.hasField("KRG") || deck.hasField("KRO") || deck.hasField("KRWR") ||
deck.hasField("KRGR") || deck.hasField("KRORW") || deck.hasField("KRORG") ||
deck.hasField("IKRW") || deck.hasField("IKRG") || deck.hasField("IKRO") || deck.hasField("IKRWR") ||
deck.hasField("IKRGR") || deck.hasField("IKRORW") || deck.hasField("IKRORG") ) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently hysteresis and relperm value scaling can not be combined.");
}
initEPSHyst(deck, grid);
}
//OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Under construction ...");
}
}
@ -165,8 +174,10 @@ namespace Opm
if (dkrds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
if (do_eps_) {
relpermEPS(s + np*i, cells[i], kr + np*i, dkrds + np*np*i);
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]]));
} else if (do_eps_) {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]]));
} else {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
}
@ -174,8 +185,10 @@ namespace Opm
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
if (do_eps_) {
relpermEPS(s + np*i, cells[i], kr + np*i);
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]]));
} else if (do_eps_) {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i, &(eps_transf_[cells[i]]));
} else {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i);
}
@ -209,12 +222,20 @@ namespace Opm
if (dpcds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i);
if (do_eps_) {
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i, &(eps_transf_[cells[i]]));
} else {
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i);
}
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i);
for (int i = 0; i < n; ++i) {
if (do_eps_) {
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i, &(eps_transf_[cells[i]]));
} else {
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i);
}
}
}
}
@ -234,16 +255,60 @@ namespace Opm
double* smax) const
{
assert(cells != 0);
const int np = phase_usage_.num_phases;
for (int i = 0; i < n; ++i) {
for (int p = 0; p < np; ++p) {
smin[np*i + p] = funcForCell(cells[i]).smin_[p];
smax[np*i + p] = funcForCell(cells[i]).smax_[p];
if (do_eps_) {
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid];
const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
for (int i = 0; i < n; ++i) {
smin[np*i + opos] = 1.0;
smax[np*i + opos] = 1.0;
if (phase_usage_.phase_used[Aqua]) {
smin[np*i + wpos] = eps_transf_[cells[i]].wat.doNotScale ? eps_transf_[cells[i]].wat.smin: funcForCell(cells[i]).smin_[wpos];
smax[np*i + wpos] = eps_transf_[cells[i]].wat.doNotScale ? eps_transf_[cells[i]].wat.smax: funcForCell(cells[i]).smax_[wpos];
smin[np*i + opos] -= smax[np*i + wpos];
smax[np*i + opos] -= smin[np*i + wpos];
}
if (phase_usage_.phase_used[Vapour]) {
smin[np*i + gpos] = eps_transf_[cells[i]].wat.doNotScale ? eps_transf_[cells[i]].gas.smin: funcForCell(cells[i]).smin_[gpos];
smax[np*i + gpos] = eps_transf_[cells[i]].wat.doNotScale ? eps_transf_[cells[i]].gas.smax: funcForCell(cells[i]).smax_[gpos];
smin[np*i + opos] -= smax[np*i + gpos];
smax[np*i + opos] -= smin[np*i + gpos];
}
if (phase_usage_.phase_used[Vapour] && phase_usage_.phase_used[Aqua]) {
smin[np*i + opos] = std::max(0.0,smin[np*i + opos]);
}
}
} else {
for (int i = 0; i < n; ++i) {
for (int p = 0; p < np; ++p) {
smin[np*i + p] = funcForCell(cells[i]).smin_[p];
smax[np*i + p] = funcForCell(cells[i]).smax_[p];
}
}
}
}
/// Update saturation state for the hysteresis tracking
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::updateSatHyst(const int n,
const int* cells,
const double* s)
{
assert(cells != 0);
const int np = phase_usage_.num_phases;
if (do_hyst_) {
// #pragma omp parallel for
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]]));
}
}
}
// Map the cell number to the correct function set.
template <class SatFuncSet>
@ -253,13 +318,147 @@ namespace Opm
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
// Initialize saturation scaling parameter
// Initialize saturation scaling parameters
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const std::string& keyword,
std::vector<double>& scaleparam)
const UnstructuredGrid& grid)
{
std::vector<double> swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr;
std::vector<double> krw, krg, kro, krwr, krgr, krorw, krorg;
// Initialize saturation scaling parameter
initEPSKey(deck, grid, std::string("SWL"), swl);
initEPSKey(deck, grid, std::string("SWU"), swu);
initEPSKey(deck, grid, std::string("SWCR"), swcr);
initEPSKey(deck, grid, std::string("SGL"), sgl);
initEPSKey(deck, grid, std::string("SGU"), sgu);
initEPSKey(deck, grid, std::string("SGCR"), sgcr);
initEPSKey(deck, grid, std::string("SOWCR"), sowcr);
initEPSKey(deck, grid, std::string("SOGCR"), sogcr);
initEPSKey(deck, grid, std::string("KRW"), krw);
initEPSKey(deck, grid, std::string("KRG"), krg);
initEPSKey(deck, grid, std::string("KRO"), kro);
initEPSKey(deck, grid, std::string("KRWR"), krwr);
initEPSKey(deck, grid, std::string("KRGR"), krgr);
initEPSKey(deck, grid, std::string("KRORW"), krorw);
initEPSKey(deck, grid, std::string("KRORG"), krorg);
eps_transf_.resize(grid.number_of_cells);
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
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 oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
for (int cell = 0; cell < grid.number_of_cells; ++cell) {
if (oilWater) {
// ### krw
initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos],
funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw);
// ### krow
initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro);
} else if (oilGas) {
// ### krg
initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos],
funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg);
// ### krog
initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro);
} else if (threephase) {
// ### krw
initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_,
funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw);
// ### krow
initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_,
funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro);
// ### krg
initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_,
funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg);
// ### krog
initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_,
funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro);
}
}
}
// Initialize hysteresis saturation scaling parameters
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(const EclipseGridParser& deck,
const UnstructuredGrid& grid)
{
std::vector<double> iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr;
std::vector<double> ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg;
// Initialize hysteresis saturation scaling parameters
initEPSKey(deck, grid, std::string("ISWL"), iswl);
initEPSKey(deck, grid, std::string("ISWU"), iswu);
initEPSKey(deck, grid, std::string("ISWCR"), iswcr);
initEPSKey(deck, grid, std::string("ISGL"), isgl);
initEPSKey(deck, grid, std::string("ISGU"), isgu);
initEPSKey(deck, grid, std::string("ISGCR"), isgcr);
initEPSKey(deck, grid, std::string("ISOWCR"), isowcr);
initEPSKey(deck, grid, std::string("ISOGCR"), isogcr);
initEPSKey(deck, grid, std::string("IKRW"), ikrw);
initEPSKey(deck, grid, std::string("IKRG"), ikrg);
initEPSKey(deck, grid, std::string("IKRO"), ikro);
initEPSKey(deck, grid, std::string("IKRWR"), ikrwr);
initEPSKey(deck, grid, std::string("IKRGR"), ikrgr);
initEPSKey(deck, grid, std::string("IKRORW"), ikrorw);
initEPSKey(deck, grid, std::string("IKRORG"), ikrorg);
eps_transf_hyst_.resize(grid.number_of_cells);
sat_hyst_.resize(grid.number_of_cells);
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
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 oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
for (int cell = 0; cell < grid.number_of_cells; ++cell) {
if (oilWater) {
// ### krw
initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos],
funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw);
// ### krow
initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro);
} else if (oilGas) {
// ### krg
initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos],
funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg);
// ### krog
initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro);
} else if (threephase) {
// ### krw
initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_,
funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw);
// ### krow
initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_,
funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro);
// ### krg
initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_,
funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg);
// ### krog
initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_,
funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro);
}
}
}
// Initialize saturation scaling parameter
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const std::string& keyword,
std::vector<double>& scaleparam)
{
const bool useAqua = phase_usage_.phase_used[Aqua];
const bool useLiquid = phase_usage_.phase_used[Liquid];
const bool useVapour = phase_usage_.phase_used[Vapour];
bool useKeyword = deck.hasField(keyword);
bool hasENPTVD = deck.hasField("ENPTVD");
bool hasENKRVD = deck.hasField("ENKRVD");
@ -269,70 +468,120 @@ namespace Opm
// Active keyword assigned default values for each cell (in case of possible box-wise assignment)
int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua];
if (keyword[0] == 'S' && (useKeyword || hasENPTVD)) {
if (keyword == std::string("SWL")) {
if (useKeyword || deck.getENPTVD().mask_[0]) {
int phase_pos_vapour = phase_usage_.phase_pos[BlackoilPhases::Vapour];
if ((keyword[0] == 'S' && (useKeyword || hasENPTVD)) || (keyword[1] == 'S' && useKeyword) ) {
if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) {
if (useAqua && (useKeyword || deck.getENPTVD().mask_[0])) {
itab = 1;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smin_[phase_pos_aqua];
}
} else if (keyword == std::string("SWCR")) {
if (useKeyword || deck.getENPTVD().mask_[1]) {
} else if (keyword == std::string("SWCR") || keyword == std::string("ISWCR") ) {
if (useAqua && (useKeyword || deck.getENPTVD().mask_[1])) {
itab = 2;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).swcr_;
}
} else if (keyword == std::string("SWU")) {
if (useKeyword || deck.getENPTVD().mask_[2]) {
} else if (keyword == std::string("SWU") || keyword == std::string("ISWU") ) {
if (useAqua && (useKeyword || deck.getENPTVD().mask_[2])) {
itab = 3;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smax_[phase_pos_aqua];
}
} else if (keyword == std::string("SOWCR")) {
if (useKeyword || deck.getENPTVD().mask_[3]) {
} else if (keyword == std::string("SGL") || keyword == std::string("ISGL") ) {
if (useVapour && (useKeyword || deck.getENPTVD().mask_[3])) {
itab = 4;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smin_[phase_pos_vapour];
}
} else if (keyword == std::string("SGCR") || keyword == std::string("ISGCR") ) {
if (useVapour && (useKeyword || deck.getENPTVD().mask_[4])) {
itab = 5;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sgcr_;
}
} else if (keyword == std::string("SGU") || keyword == std::string("ISGU") ) {
if (useVapour && (useKeyword || deck.getENPTVD().mask_[5])) {
itab = 6;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smax_[phase_pos_vapour];
}
} else if (keyword == std::string("SOWCR") || keyword == std::string("ISOWCR") ) {
if (useAqua && (useKeyword || deck.getENPTVD().mask_[6])) {
itab = 7;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sowcr_;
}
}else {
} else if (keyword == std::string("SOGCR") || keyword == std::string("ISOGCR") ) {
if (useVapour && (useKeyword || deck.getENPTVD().mask_[7])) {
itab = 8;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sogcr_;
}
} else {
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
table = deck.getENPTVD().table_;
}
} else if (keyword[0] == 'K' && (useKeyword || hasENKRVD)) {
if (keyword == std::string("KRW")) {
if (useKeyword || deck.getENKRVD().mask_[0]) {
} else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) {
if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) {
if (useAqua && (useKeyword || deck.getENKRVD().mask_[0])) {
itab = 1;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwmax_;
}
} else if (keyword == std::string("KRO")) {
if (useKeyword || deck.getENKRVD().mask_[1]) {
} else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) {
if (useVapour && (useKeyword || deck.getENKRVD().mask_[1])) {
itab = 2;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgmax_;
}
} else if (keyword == std::string("KRO") || keyword == std::string("IKRO") ) {
if (useLiquid && (useKeyword || deck.getENKRVD().mask_[2])) {
itab = 3;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).kromax_;
}
} else if (keyword == std::string("KRWR")) {
if (useKeyword || deck.getENKRVD().mask_[2]) {
itab = 3;
} else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) {
if (useAqua && (useKeyword || deck.getENKRVD().mask_[3])) {
itab = 4;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwr_;
}
} else if (keyword == std::string("KRORW")) {
if (useKeyword || deck.getENKRVD().mask_[3]) {
itab = 4;
} else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) {
if (useVapour && (useKeyword || deck.getENKRVD().mask_[4])) {
itab = 5;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgr_;
}
} else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) {
if (useAqua && (useKeyword || deck.getENKRVD().mask_[5])) {
itab = 6;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorw_;
}
} else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) {
if (useVapour && (useKeyword || deck.getENKRVD().mask_[6])) {
itab = 7;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorg_;
}
} else {
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
@ -374,171 +623,85 @@ namespace Opm
}
}
// Saturation scaling
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::relpermEPS(const double *s, const int cell, double *kr, double *dkrds) const
void SaturationPropsFromDeck<SatFuncSet>::initEPSParam(const int cell,
EPSTransforms::Transform& data,
const bool oil, // flag indicating krow/krog calculations
const double sl_tab, // minimum saturation (for krow/krog calculations this is normally zero)
const double scr_tab, // critical saturation
const double su_tab, // maximum saturation (for krow/krog calculations this is minimum water/gas saturation)
const double sxcr_tab, // second critical saturation (not used for 2pt scaling)
const double s0_tab, // threephase complementary minimum saturation (-1.0 indicates 2-phase)
const double krsr_tab, // relperm at displacing critical saturation
const double krmax_tab, // relperm at maximum saturation
const std::vector<double>& sl, // For krow/krog calculations this is not used
const std::vector<double>& scr,
const std::vector<double>& su, // For krow/krog calculations this is SWL/SGL
const std::vector<double>& sxcr,
const std::vector<double>& s0,
const std::vector<double>& krsr,
const std::vector<double>& krmax)
{
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid];
double ss[PhaseUsage::MaxNumPhases];
if (scr.empty() && su.empty() && (sxcr.empty() || !do_3pt_) && s0.empty()) {
data.doNotScale = true;
} else {
data.doNotScale = false;
data.do_3pt = do_3pt_;
double s_r;
if (s0_tab < 0.0) { // 2phase
s_r = 1.0-sxcr_tab;
if (do_3pt_) data.sr = sxcr.empty() ? s_r : 1.0-sxcr[cell];
} else { // 3phase
s_r = 1.0-sxcr_tab-s0_tab;
if (do_3pt_)data.sr = 1.0 - (sxcr.empty() ? sxcr_tab : sxcr[cell])
- (s0.empty() ? s0_tab : s0[cell]);
}
data.scr = scr.empty() ? scr_tab : scr[cell];
double s_max = su_tab;
if (oil) {
data.smin = sl_tab;
if (s0_tab < 0.0) { // 2phase
s_max = 1.0 - su_tab;
data.smax = 1.0 - (su.empty() ? su_tab : su[cell]);
} else { // 3phase
s_max = 1.0 - su_tab - s0_tab;
data.smax = 1.0 - (su.empty() ? su_tab : su[cell])
- (s0.empty() ? s0_tab : s0[cell]);
}
} else {
data.smin = sl.empty() ? sl_tab : sl[cell];
data.smax = su.empty() ? su_tab : su[cell];
}
if (do_3pt_) {
data.slope1 = (s_r-scr_tab)/(data.sr-data.scr);
data.slope2 = (s_max-s_r)/(data.smax-data.sr);
} else {
data.slope2 = data.slope1 = (s_max-scr_tab)/(data.smax-data.scr);
// Inv transform of tabulated critical displacing saturation to prepare for possible value scaling (krwr etc)
data.sr = data.scr + (s_r-scr_tab)*(data.smax-data.scr)/(s_max-scr_tab);
}
}
data.doKrMax = !krmax.empty();
data.doKrCrit = !krsr.empty();
data.doSatInterp = false;
data.krsr = krsr.empty() ? krsr_tab : krsr[cell];
data.krmax = krmax.empty() ? krmax_tab : krmax[cell];
data.krSlopeCrit = data.krsr/krsr_tab;
data.krSlopeMax = data.krmax/krmax_tab;
if (data.doKrCrit) {
if (data.sr > data.smax-1.0e-6) {
//Ignore krsr and do two-point (one might consider combining krsr and krmax linearly between scr and smax ... )
data.doKrCrit = false;
} else if (std::fabs(krmax_tab- krsr_tab) > 1.0e-6) { // interpolate kr
data.krSlopeMax = (data.krmax-data.krsr)/(krmax_tab-krsr_tab);
} else { // interpolate sat
data.doSatInterp = true;
data.krSlopeMax = (data.krmax-data.krsr)/(data.smax-data.sr);
}
}
if (do_3pt_) { // Three-point scaling
// Transforms for water saturation
if (eps_.swcr_.empty() && eps_.swu_.empty()) {
ss[wpos] = s[wpos];
} else {
double s_r = 1.0-funcForCell(cell).sowcr_;
double sr = eps_.sowcr_.empty() ? s_r : 1.0-eps_.sowcr_[cell];
if (s[wpos] <= sr) {
double sw_cr = funcForCell(cell).swcr_;
double swcr = eps_.swcr_.empty() ? sw_cr : eps_.swcr_[cell];
ss[wpos] = (s[wpos] <= swcr) ? sw_cr : sw_cr+(s[wpos]-swcr)*(s_r-sw_cr)/(sr-swcr);
} else {
double sw_max = funcForCell(cell).smax_[wpos];
double swmax = eps_.swu_.empty() ? sw_max : eps_.swu_[cell];
ss[wpos] = (s[wpos] >= swmax) ? sw_max : s_r+(s[wpos]-sr)*(sw_max-s_r)/(swmax-sr);
}
}
// Transforms for oil saturation
if (eps_.sowcr_.empty() && eps_.swl_.empty()) {
ss[opos] = s[opos];
} else {
double s_r = 1.0-funcForCell(cell).swcr_;
double sr = eps_.swcr_.empty() ? s_r : 1.0-eps_.swcr_[cell];
if (s[opos] <= sr) {
double sow_cr = funcForCell(cell).sowcr_;
double sowcr = eps_.sowcr_.empty() ? sow_cr : eps_.sowcr_[cell];
ss[opos] = (s[opos] <= sowcr) ? sow_cr : sow_cr+(s[opos]-sowcr)*(s_r-sow_cr)/(sr-sowcr);
} else {
double sow_max = funcForCell(cell).smax_[opos];
double sowmax = eps_.swl_.empty() ? sow_max : (1.0-eps_.swl_[cell]);
ss[opos] = (s[opos] >= sowmax) ? sow_max : s_r+(s[opos]-sr)*(sow_max-s_r)/(sowmax-sr);
}
}
} else { // Two-point scaling
// Transforms for water saturation
if (eps_.swcr_.empty() && eps_.swu_.empty()) {
ss[wpos] = s[wpos];
} else {
double sw_cr = funcForCell(cell).swcr_;
double swcr = eps_.swcr_.empty() ? sw_cr : eps_.swcr_[cell];
if (s[wpos] <= swcr) {
ss[wpos] = sw_cr;
} else {
double sw_max = funcForCell(cell).smax_[wpos];
double swmax = eps_.swu_.empty() ? sw_max : eps_.swu_[cell];
ss[wpos] = (s[wpos] >= swmax) ? sw_max : sw_cr + (s[wpos]-swcr)*(sw_max-sw_cr)/(swmax-swcr);
}
}
// Transforms for oil saturation
if (eps_.sowcr_.empty() && eps_.swl_.empty()) {
ss[opos] = s[opos];
} else {
double sow_cr = funcForCell(cell).sowcr_;
double socr = eps_.sowcr_.empty() ? sow_cr : eps_.sowcr_[cell];
if (s[opos] <= socr) {
ss[opos] = sow_cr;
} else {
double sow_max = funcForCell(cell).smax_[opos];
double sowmax = eps_.swl_.empty() ? sow_max : (1.0-eps_.swl_[cell]);
ss[opos] = (s[opos] >= sowmax) ? sow_max : sow_cr + (s[opos]-socr) *(sow_max-sow_cr)/(sowmax-socr);
}
}
}
// Evaluation of relperms
if (dkrds) {
OPM_THROW(std::runtime_error, "Relperm derivatives not yet available in combination with end point scaling ...");
funcForCell(cell).evalKrDeriv(ss, kr, dkrds);
} else {
// Assume: sw_cr -> krw=0 sw_max -> krw=<max water relperm>
// sow_cr -> kro=0 sow_max -> kro=<max oil relperm>
funcForCell(cell).evalKr(ss, kr);
}
// Scaling of relperms values
// - Water
if (eps_.krw_.empty() && eps_.krwr_.empty()) { // No value scaling
} else if (eps_.krwr_.empty()) { // Two-point
kr[wpos] *= (eps_.krw_[cell]/funcForCell(cell).krwmax_);
} else {
double swcr = eps_.swcr_.empty() ? funcForCell(cell).swcr_ : eps_.swcr_[cell];
double swmax = eps_.swu_.empty() ? funcForCell(cell).smax_[wpos] : eps_.swu_[cell];
double sr;
if (do_3pt_) {
sr = eps_.sowcr_.empty() ? 1.0-funcForCell(cell).sowcr_ : 1.0-eps_.sowcr_[cell];
} else {
double sw_cr = funcForCell(cell).swcr_;
double sw_max = funcForCell(cell).smax_[wpos];
double s_r = 1.0-funcForCell(cell).sowcr_;
sr = swcr + (s_r-sw_cr)*(swmax-swcr)/(sw_max-sw_cr);
}
if (s[wpos] <= swcr) {
kr[wpos] = 0.0;
} else if (sr > swmax-1.0e-6) {
if (do_3pt_) { //Ignore krw and do two-point?
kr[wpos] *= eps_.krwr_[cell]/funcForCell(cell).krwr_;
} else if (!eps_.kro_.empty()){ //Ignore krwr and do two-point
kr[wpos] *= eps_.krw_[cell]/funcForCell(cell).krwmax_;
}
} else if (s[wpos] <= sr) {
kr[wpos] *= eps_.krwr_[cell]/funcForCell(cell).krwr_;
} else if (s[wpos] <= swmax) {
double krw_max = funcForCell(cell).krwmax_;
double krw = eps_.krw_.empty() ? krw_max : eps_.krw_[cell];
double krw_r = funcForCell(cell).krwr_;
double krwr = eps_.krwr_.empty() ? krw_r : eps_.krwr_[cell];
if (std::fabs(krw_max- krw_r) > 1.0e-6) {
kr[wpos] = krwr + (kr[wpos]-krw_r)*(krw-krwr)/(krw_max-krw_r);
} else {
kr[wpos] = krwr + (krw-krwr)*(s[wpos]-sr)/(swmax-sr);
}
} else {
kr[wpos] = eps_.krw_.empty() ? funcForCell(cell).krwmax_ : eps_.krw_[cell];
}
}
// - Oil
if (eps_.kro_.empty() && eps_.krorw_.empty()) { // No value scaling
} else if (eps_.krorw_.empty()) { // Two-point scaling
kr[opos] *= (eps_.kro_[cell]/funcForCell(cell).kromax_);
} else {
double sowcr = eps_.sowcr_.empty() ? funcForCell(cell).sowcr_ : eps_.sowcr_[cell];
double sowmax = eps_.swl_.empty() ? funcForCell(cell).smax_[opos] : 1.0-eps_.swl_[cell];
double sr;
if (do_3pt_) {
sr = eps_.swcr_.empty() ? 1.0-funcForCell(cell).swcr_ : 1.0-eps_.swcr_[cell];
} else {
double sow_cr = funcForCell(cell).sowcr_;
double sow_max = funcForCell(cell).smax_[opos];
double s_r = 1.0-funcForCell(cell).swcr_;
sr = sowcr + (s_r-sow_cr)*(sowmax-sowcr)/(sow_max-sow_cr);
}
if (s[opos] <= sowcr) {
kr[opos] = 0.0;
} else if (sr > sowmax-1.0e-6) {
if (do_3pt_) { //Ignore kro and do two-point?
kr[opos] *= eps_.krorw_[cell]/funcForCell(cell).krorw_;
} else if (!eps_.kro_.empty()){ //Ignore krowr and do two-point
kr[opos] *= eps_.kro_[cell]/funcForCell(cell).kromax_;
}
} else if (s[opos] <= sr) {
kr[opos] *= eps_.krorw_[cell]/funcForCell(cell).krorw_;
} else if (s[opos] <= sowmax) {
double kro_max = funcForCell(cell).kromax_;
double kro = eps_.kro_.empty() ? kro_max : eps_.kro_[cell];
double kro_rw = funcForCell(cell).krorw_;
double krorw = eps_.krorw_[cell];
if (std::fabs(kro_max- kro_rw) > 1.0e-6) {
kr[opos] = krorw + (kr[opos]- kro_rw)*(kro-krorw)/(kro_max- kro_rw);
} else {
kr[opos] = krorw + (kro-krorw)*(s[opos]-sr)/(sowmax-sr);
}
} else {
kr[opos] = eps_.kro_.empty() ? funcForCell(cell).kromax_ : eps_.kro_[cell];
}
}
}
} // namespace Opm

3
tests/testBlackoilState1.DATA
View File

@ -13,3 +13,6 @@ DZV
ACTNUM
1 998*2 3 /
DEPTHZ
121*2000 /

3
tests/testBlackoilState2.DATA
View File

@ -9,3 +9,6 @@ DYV
DZV
10*10 /
DEPTHZ
110*2000 /

2
tests/test_parser.cpp
View File

@ -43,7 +43,7 @@ BOOST_AUTO_TEST_CASE(CreateParser)
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile( filename1 );
BOOST_CHECK_EQUAL( 5U , deck->size() );
BOOST_CHECK_EQUAL( 6U , deck->size() );
Opm::DeckItemConstPtr actnum = deck->getKeyword("ACTNUM")->getRecord(0)->getItem(0);
const std::vector<int>& actnum_data = actnum->getIntData();

4
tests/wells_manager_data.data
View File

@ -16,6 +16,10 @@ DYV
DZV
10.0 20.0 30.0 10.0 5.0 /
DEPTHZ
121*2000
/
SCHEDULE
WELSPECS

6
tests/wells_manager_data_expanded.data
View File

@ -16,6 +16,12 @@ DYV
DZV
10.0 20.0 30.0 10.0 5.0 /
-- The DEPTHZ keyword is only here to satisfy the old parser; content might
-- completely bogus.
DEPTHZ
121*2000 /
SCHEDULE
WELSPECS