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

Merge pull request #652 from andlaus/table_improvements

Browse Source 
adapt the the table related API changes of opm-parser
This commit is contained in:
Joakim Hove 2014-09-19 15:27:29 +02:00
commit cb3b1c18fe
24 changed files with 501 additions and 493 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/sim_2p_comp_reorder.cpp
View File

@ -108,7 +108,7 @@ try
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
rock_comp.reset(new RockCompressibility(deck));
rock_comp.reset(new RockCompressibility(deck, eclipseState));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).

2
examples/sim_2p_incomp.cpp
View File

@ -121,7 +121,7 @@ try
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
rock_comp.reset(new RockCompressibility(deck));
rock_comp.reset(new RockCompressibility(deck, eclipseState));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).

2
examples/wells_example.cpp
View File

@ -44,7 +44,7 @@ try
// Define rock and fluid properties
IncompPropertiesFromDeck incomp_properties(deck, eclipseState, *grid.c_grid());
RockCompressibility rock_comp(deck);
RockCompressibility rock_comp(deck, eclipseState);
// Finally handle the wells
WellsManager wells(eclipseState , 0 , *grid.c_grid(), incomp_properties.permeability());

4
opm/core/props/BlackoilPropertiesFromDeck_impl.hpp
View File

@ -53,7 +53,7 @@ namespace Opm
if (init_rock){
rock_.init(eclState, number_of_cells, global_cell, cart_dims);
}
pvt_.init(deck, /*numSamples=*/0);
pvt_.init(deck, eclState, /*numSamples=*/0);
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
@ -86,7 +86,7 @@ namespace Opm
}
const int pvt_samples = param.getDefault("pvt_tab_size", -1);
pvt_.init(deck, pvt_samples);
pvt_.init(deck, eclState, pvt_samples);
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", -1);

31
opm/core/props/pvt/BlackoilPvtProperties.cpp
View File

@ -31,9 +31,8 @@
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
#include <opm/parser/eclipse/Utility/PvcdoTable.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
namespace Opm
{
@ -43,6 +42,7 @@ namespace Opm
}
void BlackoilPvtProperties::init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int numSamples)
{
phase_usage_ = phaseUsageFromDeck(deck);
@ -83,19 +83,20 @@ namespace Opm
if (phase_usage_.phase_used[Liquid]) {
// for oil, we support the "PVDO", "PVTO" and "PVCDO"
// keywords...
if (deck->hasKeyword("PVDO")) {
Opm::DeckKeywordConstPtr pvdoKeyword = deck->getKeyword("PVDO");
const auto &pvdoTables = eclipseState->getPvdoTables();
const auto &pvtoTables = eclipseState->getPvtoTables();
if (pvdoTables.size() > 0) {
if (numSamples > 0) {
auto splinePvt = std::shared_ptr<PvtDeadSpline>(new PvtDeadSpline);
splinePvt->initFromOil(pvdoKeyword, numSamples);
splinePvt->initFromOil(pvdoTables, numSamples);
props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
} else {
auto deadPvt = std::shared_ptr<PvtDead>(new PvtDead);
deadPvt->initFromOil(pvdoKeyword);
deadPvt->initFromOil(pvdoTables);
props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
}
} else if (deck->hasKeyword("PVTO")) {
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO")));
} else if (pvtoTables.size() > 0) {
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(pvtoTables));
} else if (deck->hasKeyword("PVCDO")) {
std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
pvcdo->initFromOil(deck->getKeyword("PVCDO"));
@ -108,22 +109,22 @@ namespace Opm
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
// gas can be specified using the "PVDG" or "PVTG" keywords...
if (deck->hasKeyword("PVDG")) {
Opm::DeckKeywordConstPtr pvdgKeyword = deck->getKeyword("PVDG");
const auto &pvdgTables = eclipseState->getPvdgTables();
const auto &pvtgTables = eclipseState->getPvtgTables();
if (pvdgTables.size() > 0) {
if (numSamples > 0) {
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromGas(pvdgKeyword, numSamples);
splinePvt->initFromGas(pvdgTables, numSamples);
props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
} else {
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromGas(pvdgKeyword);
deadPvt->initFromGas(pvdgTables);
props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
}
} else if (deck->hasKeyword("PVTG")) {
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG")));
} else if (pvtgTables.size() > 0) {
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(pvtgTables));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
}

3
opm/core/props/pvt/BlackoilPvtProperties.hpp
View File

@ -24,7 +24,7 @@
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <string>
#include <memory>
@ -54,6 +54,7 @@ namespace Opm
///
/// \param deck An input deck from the opm-parser module.
void init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int samples);
/// \return Object describing the active phases.

3
opm/core/props/pvt/PvtConstCompr.hpp
View File

@ -23,9 +23,6 @@
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
#include <opm/parser/eclipse/Utility/PvcdoTable.hpp>
#include <vector>
#include <algorithm>

13
opm/core/props/pvt/PvtDead.cpp
View File

@ -34,16 +34,17 @@ namespace Opm
// Member functions
//-------------------------------------------------------------------------
/// Constructor
void PvtDead::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword)
void PvtDead::initFromOil(const std::vector<Opm::PvdoTable>& pvdoTables)
{
int numRegions = Opm::PvdoTable::numTables(pvdoKeyword);
int numRegions = pvdoTables.size();
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx);
const Opm::PvdoTable& pvdoTable = pvdoTables[regionIdx];
// Copy data
const std::vector<double>& press = pvdoTable.getPressureColumn();
const std::vector<double>& b = pvdoTable.getFormationFactorColumn();
@ -60,16 +61,16 @@ namespace Opm
}
void PvtDead::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword)
void PvtDead::initFromGas(const std::vector<Opm::PvdgTable>& pvdgTables)
{
int numRegions = Opm::PvdgTable::numTables(pvdgKeyword);
int numRegions = pvdgTables.size();
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx);
const Opm::PvdgTable& pvdgTable = pvdgTables[regionIdx];
// Copy data
const std::vector<double>& press = pvdgTable.getPressureColumn();

8
opm/core/props/pvt/PvtDead.hpp
View File

@ -20,12 +20,10 @@
#ifndef OPM_PVTDEAD_HEADER_INCLUDED
#define OPM_PVTDEAD_HEADER_INCLUDED
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
#include <opm/parser/eclipse/Utility/PvdgTable.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
@ -48,8 +46,8 @@ namespace Opm
public:
PvtDead() {};
void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword);
void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword);
void initFromOil(const std::vector<Opm::PvdoTable>& pvdoTables);
void initFromGas(const std::vector<Opm::PvdgTable>& pvdgTables);
virtual ~PvtDead();
/// Viscosity as a function of p and z.

13
opm/core/props/pvt/PvtDeadSpline.cpp
View File

@ -21,7 +21,6 @@
#include <opm/core/props/pvt/PvtDeadSpline.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/parser/eclipse/Utility/SingleRecordTable.hpp>
#include <algorithm>
@ -39,17 +38,17 @@ namespace Opm
PvtDeadSpline::PvtDeadSpline()
{}
void PvtDeadSpline::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
void PvtDeadSpline::initFromOil(const std::vector<Opm::PvdoTable>& pvdoTables,
int numSamples)
{
int numRegions = Opm::PvdoTable::numTables(pvdoKeyword);
int numRegions = pvdoTables.size();
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx);
const Opm::PvdoTable& pvdoTable = pvdoTables[regionIdx];
int numRows = pvdoTable.numRows();
@ -68,17 +67,17 @@ namespace Opm
}
}
void PvtDeadSpline::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
void PvtDeadSpline::initFromGas(const std::vector<Opm::PvdgTable>& pvdgTables,
int numSamples)
{
int numRegions = Opm::PvdgTable::numTables(pvdgKeyword);
int numRegions = pvdgTables.size();
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx);
const Opm::PvdgTable& pvdgTable = pvdgTables[regionIdx];
int numRows = pvdgTable.numRows();

7
opm/core/props/pvt/PvtDeadSpline.hpp
View File

@ -23,8 +23,7 @@
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
#include <opm/parser/eclipse/Utility/PvdgTable.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
@ -43,9 +42,9 @@ namespace Opm
public:
PvtDeadSpline();
void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
void initFromOil(const std::vector<Opm::PvdoTable>& pvdoTables,
int numSamples);
void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
void initFromGas(const std::vector<Opm::PvdgTable>& pvdgTables,
int numSamples);
virtual ~PvtDeadSpline();

6
opm/core/props/pvt/PvtLiveGas.cpp
View File

@ -46,14 +46,14 @@ namespace Opm
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
PvtLiveGas::PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword)
PvtLiveGas::PvtLiveGas(const std::vector<Opm::PvtgTable>& pvtgTables)
{
int numTables = Opm::PvtgTable::numTables(pvtgKeyword);
int numTables = pvtgTables.size();
saturated_gas_table_.resize(numTables);
undersat_gas_tables_.resize(numTables);
for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) {
Opm::PvtgTable pvtgTable(pvtgKeyword, pvtTableIdx);
const Opm::PvtgTable& pvtgTable = pvtgTables[pvtTableIdx];
// GAS, PVTG
saturated_gas_table_[pvtTableIdx].resize(4);

4
opm/core/props/pvt/PvtLiveGas.hpp
View File

@ -22,7 +22,7 @@
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/parser/eclipse/Utility/PvtgTable.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
@ -38,7 +38,7 @@ namespace Opm
class PvtLiveGas : public PvtInterface
{
public:
PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword);
PvtLiveGas(const std::vector<Opm::PvtgTable>& pvtgTables);
virtual ~PvtLiveGas();
/// Viscosity as a function of p and z.

6
opm/core/props/pvt/PvtLiveOil.cpp
View File

@ -34,14 +34,14 @@ namespace Opm
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
PvtLiveOil::PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword)
PvtLiveOil::PvtLiveOil(const std::vector<Opm::PvtoTable>& pvtoTables)
{
int numTables = Opm::PvtoTable::numTables(pvtoKeyword);
int numTables = pvtoTables.size();
saturated_oil_table_.resize(numTables);
undersat_oil_tables_.resize(numTables);
for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) {
Opm::PvtoTable pvtoTable(pvtoKeyword, pvtTableIdx);
Opm::PvtoTable pvtoTable = pvtoTables[pvtTableIdx];
const auto saturatedPvto = pvtoTable.getOuterTable();

4
opm/core/props/pvt/PvtLiveOil.hpp
View File

@ -23,7 +23,7 @@
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
@ -39,7 +39,7 @@ namespace Opm
class PvtLiveOil : public PvtInterface
{
public:
PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword);
PvtLiveOil(const std::vector<Opm::PvtoTable>& pvtoTables);
virtual ~PvtLiveOil();
/// Viscosity as a function of p and z.

45
opm/core/props/rock/RockCompressibility.cpp
View File

@ -24,8 +24,6 @@
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/Utility/RocktabTable.hpp>
#include <opm/parser/eclipse/Utility/RockTable.hpp>
#include <iostream>
@ -40,39 +38,30 @@ namespace Opm
rock_comp_ = param.getDefault("rock_compressibility", 0.0)/unit::barsa;
}
RockCompressibility::RockCompressibility(Opm::DeckConstPtr deck)
RockCompressibility::RockCompressibility(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState)
: pref_(0.0),
rock_comp_(0.0)
{
if (deck->hasKeyword("ROCKTAB")) {
Opm::DeckKeywordConstPtr rtKeyword = deck->getKeyword("ROCKTAB");
if (rtKeyword->size() != 1)
const auto& rocktabTables = eclipseState->getRocktabTables();
if (rocktabTables.size() > 0) {
if (rocktabTables.size() != 1)
OPM_THROW(std::runtime_error, "Can only handle a single region in ROCKTAB.");
// the number of colums of the "ROCKTAB" keyword
// depends on the presence of the "RKTRMDIR"
// keyword. Messy stuff...
bool isDirectional = deck->hasKeyword("RKTRMDIR");
if (isDirectional)
{
// well, okay. we don't support non-isotropic
// transmissibility multipliers yet
OPM_THROW(std::runtime_error, "Support for non-isotropic "
"transmissibility multipliers is not implemented yet.");
};
Opm::RocktabTable rocktabTable(rtKeyword, isDirectional);
p_ = rocktabTable.getPressureColumn();
poromult_ = rocktabTable.getPoreVolumeMultiplierColumn();
transmult_ = rocktabTable.getTransmissibilityMultiplierColumn();
p_ = rocktabTables[0].getPressureColumn();
poromult_ = rocktabTables[0].getPoreVolumeMultiplierColumn();
transmult_ = rocktabTables[0].getTransmissibilityMultiplierColumn();
} else if (deck->hasKeyword("ROCK")) {
Opm::RockTable rockTable(deck->getKeyword("ROCK"));
if (rockTable.numRows() != 1)
OPM_THROW(std::runtime_error, "Can only handle a single region in ROCK.");
Opm::DeckKeywordConstPtr rockKeyword = deck->getKeyword("ROCK");
if (rockKeyword->size() != 1) {
// here it would be better not to use std::cout directly but to add the
// warning to some "warning list"...
std::cout << "Can only handle a single region in ROCK ("<<rockKeyword->size()<<" regions specified)."
<< " Ignoring all except for the first.\n";
}
pref_ = rockTable.getPressureColumn()[0];
rock_comp_ = rockTable.getCompressibilityColumn()[0];
pref_ = rockKeyword->getRecord(0)->getItem("PREF")->getSIDouble(0);
rock_comp_ = rockKeyword->getRecord(0)->getItem("COMPRESSIBILITY")->getSIDouble(0);
} else {
std::cout << "**** warning: no rock compressibility data found in deck (ROCK or ROCKTAB)." << std::endl;
}

4
opm/core/props/rock/RockCompressibility.hpp
View File

@ -21,6 +21,7 @@
#define OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
@ -34,7 +35,8 @@ namespace Opm
public:
/// Construct from input deck.
/// Looks for the keywords ROCK and ROCKTAB.
RockCompressibility(Opm::DeckConstPtr deck);
RockCompressibility(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState);
/// Construct from parameters.
/// Accepts the following parameters (with defaults).

640
opm/core/props/satfunc/SatFuncBase.hpp
View File

@ -1,321 +1,319 @@
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef SATFUNCBASE_HPP
#define SATFUNCBASE_HPP
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Utility/SwofTable.hpp>
#include <opm/parser/eclipse/Utility/SgofTable.hpp>
#include <vector>
namespace Opm
{
// Transforms for saturation table scaling
struct EPSTransforms {
struct Transform {
bool doNotScale;
bool do_3pt;
double smin;
double scr;
double sr;
double smax;
double slope1;
double slope2;
double scaleSat(double ss, double s_r, double s_cr, double s_max) const;
double scaleSatInv(double s, double s_r, double s_cr, double s_max) const;
double scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const; // Returns scaleSat'(s)
double scaleSatPc(double s, double s_min, double s_max) const;
double scaleSatDerivPc(double s, double s_min, double s_max) const; // Returns scaleSatPc'(s)
bool doKrMax;
bool doKrCrit;
bool doSatInterp;
double krsr;
double krmax;
double krSlopeMax;
double krSlopeCrit;
double scaleKr(double s, double kr, double krsr_tab) const;
double scaleKrDeriv(double s, double krDeriv) const; // Returns scaleKr'(kr(scaleSat(s)))*kr'((scaleSat(s))
double pcFactor; // Scaling factor for capillary pressure.
void printMe(std::ostream & out);
};
Transform wat;
Transform watoil;
Transform gas;
Transform gasoil;
};
// Hysteresis
struct SatHyst {
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(Opm::DeckConstPtr 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.
double pcwmax_; // Max oil-water capillary pressure.
double pcgmax_; // Max gas-oil capillary pressure.
protected:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
TableType krw_;
TableType krow_;
TableType pcow_;
TableType krg_;
TableType krog_;
TableType pcog_;
double krocw_; // = krow_(s_wc)
private:
void extendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
void initializeTableType(TableType& table,
const std::vector<double>& arg,
const std::vector<double>& value,
const int samples);
};
template <class TableType>
void SatFuncBase<TableType>::init(Opm::DeckConstPtr 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]) {
Opm::SwofTable swof(deck->getKeyword("SWOF"), table_num);
const std::vector<double>& sw = swof.getSwColumn();
const std::vector<double>& krw = swof.getKrwColumn();
const std::vector<double>& krow = swof.getKrowColumn();
const std::vector<double>& pcow = swof.getPcowColumn();
if (krw.front() != 0.0 || krow.back() != 0.0) {
OPM_THROW(std::runtime_error, "Error SWOF data - non-zero krw(swco) and/or krow(1-sor)");
}
// Extend the tables with constant values such that the
// derivatives at the endpoints are zero
int n = sw.size();
std::vector<double> sw_ex(n+2);
std::vector<double> krw_ex(n+2);
std::vector<double> krow_ex(n+2);
std::vector<double> pcow_ex(n+2);
extendTable(sw,sw_ex,1);
extendTable(krw,krw_ex,0);
extendTable(krow,krow_ex,0);
extendTable(pcow,pcow_ex,0);
initializeTableType(krw_,sw_ex, krw_ex, samples);
initializeTableType(krow_,sw_ex, krow_ex, samples);
initializeTableType(pcow_,sw_ex, pcow_ex, samples);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
krwmax_ = krw.back();
kromax_ = krow.front();
swcr_ = swmax;
sowcr_ = 1.0 - swco;
krwr_ = krw.back();
krorw_ = krow.front();
for (std::vector<double>::size_type i=1; i<sw.size(); ++i) {
if (krw[i]> 0.0) {
swcr_ = sw[i-1];
krorw_ = krow[i-1];
break;
}
}
for (std::vector<double>::size_type i=sw.size()-1; i>=1; --i) {
if (krow[i-1]> 0.0) {
sowcr_ = 1.0 - sw[i];
krwr_ = krw[i];
break;
}
}
pcwmax_ = pcow.front();
}
if (phase_usage.phase_used[Vapour]) {
Opm::SgofTable sgof(deck->getKeyword("SGOF"), table_num);
const std::vector<double>& sg = sgof.getSgColumn();
const std::vector<double>& krg = sgof.getKrgColumn();
const std::vector<double>& krog = sgof.getKrogColumn();
const std::vector<double>& pcog = sgof.getPcogColumn();
// Extend the tables with constant values such that the
// derivatives at the endpoints are zero
int n = sg.size();
std::vector<double> sg_ex(n+2);
std::vector<double> krg_ex(n+2);
std::vector<double> krog_ex(n+2);
std::vector<double> pcog_ex(n+2);
extendTable(sg,sg_ex,1);
extendTable(krg,krg_ex,0);
extendTable(krog,krog_ex,0);
extendTable(pcog,pcog_ex,0);
initializeTableType(krg_,sg_ex, krg_ex, samples);
initializeTableType(krog_,sg_ex, krog_ex, samples);
initializeTableType(pcog_,sg_ex, pcog_ex, samples);
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
smin_[phase_usage.phase_pos[Vapour]] = sg.front();
krgmax_ = krg.back();
sgcr_ = sg.front();
sogcr_ = 1.0 - sg.back();
krgr_ = krg.back();
krorg_ = krg.front();
for (std::vector<double>::size_type i=1; i<sg.size(); ++i) {
if (krg[i]> 0.0) {
sgcr_ = sg[i-1];
krorg_ = krog[i-1];
break;
}
}
for (std::vector<double>::size_type i=sg.size()-1; i>=1; --i) {
if (krog[i-1]> 0.0) {
sogcr_ = 1.0 - sg[i];
krgr_ = krg[i];
break;
}
}
pcgmax_ = pcog.back();
}
if (phase_usage.phase_used[Vapour] && phase_usage.phase_used[Aqua]) {
sowcr_ -= smin_[phase_usage.phase_pos[Vapour]];
sogcr_ -= smin_[phase_usage.phase_pos[Aqua]];
smin_[phase_usage.phase_pos[Liquid]] = 0.0;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]]
- smin_[phase_usage.phase_pos[Vapour]]; // First entry in SGOF-table supposed to be zero anyway ...
} else if (phase_usage.phase_used[Aqua]) {
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Aqua]];
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]];
} else if (phase_usage.phase_used[Vapour]) {
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Vapour]];
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Vapour]];
}
}
template <class TableType>
void SatFuncBase<TableType>::updateSatHyst(const double* s,
const EPSTransforms* epst,
const EPSTransforms* epst_hyst,
SatHyst* sat_hyst) const
{
if (phase_usage.phase_used[Aqua] && phase_usage.phase_used[Vapour]) { //Water/Oil/Gas
int opos = phase_usage.phase_pos[BlackoilPhases::Liquid];
int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour];
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
if (s[opos] > sat_hyst->sow_hyst)
{
sat_hyst->sow_hyst = s[opos];
double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]);
double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]);
sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst;
}
if (s[gpos] > sat_hyst->sg_hyst)
{
sat_hyst->sg_hyst = s[gpos];
double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]);
double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]);
sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst;
}
} else if (phase_usage.phase_used[Aqua]) { //Water/oil
int opos = phase_usage.phase_pos[BlackoilPhases::Liquid];
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
if (s[opos] > sat_hyst->sow_hyst)
{
sat_hyst->sow_hyst = s[opos];
double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]);
double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]);
sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst;
}
} else if (phase_usage.phase_used[Vapour]) {//Gas/Oil
int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour];
if (s[gpos] > sat_hyst->sg_hyst)
{
sat_hyst->sg_hyst = s[gpos];
double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]);
double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]);
sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst;
}
}
}
template <class TableType>
void SatFuncBase<TableType>::extendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const
{
int n = xv.size();
xv_ex[0] = xv[0]-pm;
xv_ex[n+1] = xv[n-1]+pm;
for (int i=0; i<n; i++)
{
xv_ex[i+1] = xv[i];
}
}
} // namespace Opm
#endif // SATFUNCBASE_HPP
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef SATFUNCBASE_HPP
#define SATFUNCBASE_HPP
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
namespace Opm
{
// Transforms for saturation table scaling
struct EPSTransforms {
struct Transform {
bool doNotScale;
bool do_3pt;
double smin;
double scr;
double sr;
double smax;
double slope1;
double slope2;
double scaleSat(double ss, double s_r, double s_cr, double s_max) const;
double scaleSatInv(double s, double s_r, double s_cr, double s_max) const;
double scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const; // Returns scaleSat'(s)
double scaleSatPc(double s, double s_min, double s_max) const;
double scaleSatDerivPc(double s, double s_min, double s_max) const; // Returns scaleSatPc'(s)
bool doKrMax;
bool doKrCrit;
bool doSatInterp;
double krsr;
double krmax;
double krSlopeMax;
double krSlopeCrit;
double scaleKr(double s, double kr, double krsr_tab) const;
double scaleKrDeriv(double s, double krDeriv) const; // Returns scaleKr'(kr(scaleSat(s)))*kr'((scaleSat(s))
double pcFactor; // Scaling factor for capillary pressure.
void printMe(std::ostream & out);
};
Transform wat;
Transform watoil;
Transform gas;
Transform gasoil;
};
// Hysteresis
struct SatHyst {
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(Opm::EclipseStateConstPtr eclipseState,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void updateSatHyst(const double* s,
const EPSTransforms* epst,
const EPSTransforms* epst_hyst,
SatHyst* sat_hyst) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double krgmax_; // Max gas relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double sgcr_; // Critical gas saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double krgr_; // Gas relperm at critical oil-in-gas saturation.
double sowcr_; // Critical oil-in-water saturation.
double sogcr_; // Critical oil-in-gas-and-connate-water saturation.
double krorw_; // Oil relperm at critical water saturation.
double krorg_; // Oil relperm at critical gas saturation.
double pcwmax_; // Max oil-water capillary pressure.
double pcgmax_; // Max gas-oil capillary pressure.
protected:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
TableType krw_;
TableType krow_;
TableType pcow_;
TableType krg_;
TableType krog_;
TableType pcog_;
double krocw_; // = krow_(s_wc)
private:
void extendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
void initializeTableType(TableType& table,
const std::vector<double>& arg,
const std::vector<double>& value,
const int samples);
};
template <class TableType>
void SatFuncBase<TableType>::init(Opm::EclipseStateConstPtr eclipseState,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const Opm::SwofTable& swof(eclipseState->getSwofTables()[table_num]);
const std::vector<double>& sw = swof.getSwColumn();
const std::vector<double>& krw = swof.getKrwColumn();
const std::vector<double>& krow = swof.getKrowColumn();
const std::vector<double>& pcow = swof.getPcowColumn();
if (krw.front() != 0.0 || krow.back() != 0.0) {
OPM_THROW(std::runtime_error, "Error SWOF data - non-zero krw(swco) and/or krow(1-sor)");
}
// Extend the tables with constant values such that the
// derivatives at the endpoints are zero
int n = sw.size();
std::vector<double> sw_ex(n+2);
std::vector<double> krw_ex(n+2);
std::vector<double> krow_ex(n+2);
std::vector<double> pcow_ex(n+2);
extendTable(sw,sw_ex,1);
extendTable(krw,krw_ex,0);
extendTable(krow,krow_ex,0);
extendTable(pcow,pcow_ex,0);
initializeTableType(krw_,sw_ex, krw_ex, samples);
initializeTableType(krow_,sw_ex, krow_ex, samples);
initializeTableType(pcow_,sw_ex, pcow_ex, samples);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
krwmax_ = krw.back();
kromax_ = krow.front();
swcr_ = swmax;
sowcr_ = 1.0 - swco;
krwr_ = krw.back();
krorw_ = krow.front();
for (std::vector<double>::size_type i=1; i<sw.size(); ++i) {
if (krw[i]> 0.0) {
swcr_ = sw[i-1];
krorw_ = krow[i-1];
break;
}
}
for (std::vector<double>::size_type i=sw.size()-1; i>=1; --i) {
if (krow[i-1]> 0.0) {
sowcr_ = 1.0 - sw[i];
krwr_ = krw[i];
break;
}
}
pcwmax_ = pcow.front();
}
if (phase_usage.phase_used[Vapour]) {
const Opm::SgofTable& sgof = eclipseState->getSgofTables()[table_num];
const std::vector<double>& sg = sgof.getSgColumn();
const std::vector<double>& krg = sgof.getKrgColumn();
const std::vector<double>& krog = sgof.getKrogColumn();
const std::vector<double>& pcog = sgof.getPcogColumn();
// Extend the tables with constant values such that the
// derivatives at the endpoints are zero
int n = sg.size();
std::vector<double> sg_ex(n+2);
std::vector<double> krg_ex(n+2);
std::vector<double> krog_ex(n+2);
std::vector<double> pcog_ex(n+2);
extendTable(sg,sg_ex,1);
extendTable(krg,krg_ex,0);
extendTable(krog,krog_ex,0);
extendTable(pcog,pcog_ex,0);
initializeTableType(krg_,sg_ex, krg_ex, samples);
initializeTableType(krog_,sg_ex, krog_ex, samples);
initializeTableType(pcog_,sg_ex, pcog_ex, samples);
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
smin_[phase_usage.phase_pos[Vapour]] = sg.front();
krgmax_ = krg.back();
sgcr_ = sg.front();
sogcr_ = 1.0 - sg.back();
krgr_ = krg.back();
krorg_ = krg.front();
for (std::vector<double>::size_type i=1; i<sg.size(); ++i) {
if (krg[i]> 0.0) {
sgcr_ = sg[i-1];
krorg_ = krog[i-1];
break;
}
}
for (std::vector<double>::size_type i=sg.size()-1; i>=1; --i) {
if (krog[i-1]> 0.0) {
sogcr_ = 1.0 - sg[i];
krgr_ = krg[i];
break;
}
}
pcgmax_ = pcog.back();
}
if (phase_usage.phase_used[Vapour] && phase_usage.phase_used[Aqua]) {
sowcr_ -= smin_[phase_usage.phase_pos[Vapour]];
sogcr_ -= smin_[phase_usage.phase_pos[Aqua]];
smin_[phase_usage.phase_pos[Liquid]] = 0.0;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]]
- smin_[phase_usage.phase_pos[Vapour]]; // First entry in SGOF-table supposed to be zero anyway ...
} else if (phase_usage.phase_used[Aqua]) {
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Aqua]];
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]];
} else if (phase_usage.phase_used[Vapour]) {
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Vapour]];
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Vapour]];
}
}
template <class TableType>
void SatFuncBase<TableType>::updateSatHyst(const double* s,
const EPSTransforms* epst,
const EPSTransforms* epst_hyst,
SatHyst* sat_hyst) const
{
if (phase_usage.phase_used[Aqua] && phase_usage.phase_used[Vapour]) { //Water/Oil/Gas
int opos = phase_usage.phase_pos[BlackoilPhases::Liquid];
int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour];
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
if (s[opos] > sat_hyst->sow_hyst)
{
sat_hyst->sow_hyst = s[opos];
double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]);
double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]);
sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst;
}
if (s[gpos] > sat_hyst->sg_hyst)
{
sat_hyst->sg_hyst = s[gpos];
double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]);
double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]);
sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst;
}
} else if (phase_usage.phase_used[Aqua]) { //Water/oil
int opos = phase_usage.phase_pos[BlackoilPhases::Liquid];
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
if (s[opos] > sat_hyst->sow_hyst)
{
sat_hyst->sow_hyst = s[opos];
double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]);
double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]);
sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst;
}
} else if (phase_usage.phase_used[Vapour]) {//Gas/Oil
int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour];
if (s[gpos] > sat_hyst->sg_hyst)
{
sat_hyst->sg_hyst = s[gpos];
double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]);
double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]);
sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst;
}
}
}
template <class TableType>
void SatFuncBase<TableType>::extendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const
{
int n = xv.size();
xv_ex[0] = xv[0]-pm;
xv_ex[n+1] = xv[n-1]+pm;
for (int i=0; i<n; i++)
{
xv_ex[i+1] = xv[i];
}
}
} // namespace Opm
#endif // SATFUNCBASE_HPP

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

@ -161,18 +161,21 @@ namespace Opm
const Funcs& funcForCell(const int cell) const;
template<class T>
void initEPS(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSHyst(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSKey(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,

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

@ -29,8 +29,6 @@
#include <opm/parser/eclipse/Utility/EndscaleWrapper.hpp>
#include <opm/parser/eclipse/Utility/ScalecrsWrapper.hpp>
#include <opm/parser/eclipse/Utility/EnptvdTable.hpp>
#include <opm/parser/eclipse/Utility/EnkrvdTable.hpp>
#include <iostream>
@ -63,7 +61,7 @@ namespace Opm
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr /*eclState*/,
Opm::EclipseStateConstPtr eclState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
@ -130,7 +128,7 @@ namespace Opm
// Initialize tables.
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_, samples);
satfuncset_[table].init(eclState, table, phase_usage_, samples);
}
// Check EHYSTR status
@ -192,7 +190,7 @@ namespace Opm
// TODO: ENPTVD/ENKRVD: Too few tables gives a cryptical message from parser,
// superfluous tables are ignored by the parser without any warning ...
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids,
initEPS(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
dimensions);
if (do_hyst_) {
@ -233,7 +231,7 @@ namespace Opm
// to be a scaled version of the drainage curve (confer Norne model).
}
initEPSHyst(deck, number_of_cells, global_cell, begin_cell_centroids,
initEPSHyst(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
dimensions);
}
}
@ -460,6 +458,7 @@ namespace Opm
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
@ -470,39 +469,39 @@ namespace Opm
std::vector<double> pcw, pcg;
const std::vector<double> dummy;
// Initialize saturation scaling parameter
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWL"), swl);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWU"), swu);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWCR"), swcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGL"), sgl);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGU"), sgu);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGCR"), sgcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SOWCR"), sowcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SOGCR"), sogcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRW"), krw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRG"), krg);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRO"), kro);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRWR"), krwr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRGR"), krgr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRORW"), krorw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRORG"), krorg);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("PCW"), pcw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("PCG"), pcg);
eps_transf_.resize(number_of_cells);
@ -613,6 +612,7 @@ namespace Opm
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
@ -623,39 +623,39 @@ namespace Opm
std::vector<double> ipcw, ipcg;
const std::vector<double> dummy;
// Initialize hysteresis saturation scaling parameters
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWL"), iswl);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWU"), iswu);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWCR"), iswcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGL"), isgl);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGU"), isgu);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGCR"), isgcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISOWCR"), isowcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISOGCR"), isogcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRW"), ikrw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRG"), ikrg);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRO"), ikro);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRWR"), ikrwr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRGR"), ikrgr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRORW"), ikrorw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRORG"), ikrorg);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IPCW"), ipcw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IPCG"), ipcg);
eps_transf_hyst_.resize(number_of_cells);
@ -767,6 +767,7 @@ namespace Opm
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
@ -849,14 +850,15 @@ namespace Opm
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
int num_tables = deck->getKeyword("ENPTVD")->size();
const auto& enptvdTables = eclipseState->getEnptvdTables();
int num_tables = enptvdTables.size();
param_col.resize(num_tables);
depth_col.resize(num_tables);
col_names.resize(9);
for (int table_num=0; table_num<num_tables; ++table_num) {
Opm::SingleRecordTable enptvd(deck->getKeyword("ENPTVD"), col_names, table_num);
depth_col[table_num] = enptvd.getColumn(0); // depth
param_col[table_num] = enptvd.getColumn(itab); // itab=[1-8]: swl swcr swu sgl sgcr sgu sowcr sogcr
const auto& enptvdTable = enptvdTables[table_num];
depth_col[table_num] = enptvdTable.getDepthColumn();
param_col[table_num] = enptvdTable.getColumn(itab); // itab=[1-8]: swl swcr swu sgl sgcr sgu sowcr sogcr
}
}
} else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) {
@ -913,14 +915,15 @@ namespace Opm
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
int num_tables = deck->getKeyword("ENKRVD")->size();
const auto& enkrvdTables = eclipseState->getEnkrvdTables();
int num_tables = enkrvdTables.size();
param_col.resize(num_tables);
depth_col.resize(num_tables);
col_names.resize(8);
for (int table_num=0; table_num<num_tables; ++table_num) {
Opm::SingleRecordTable enkrvd(deck->getKeyword("ENKRVD"), col_names, table_num);
depth_col[table_num] = enkrvd.getColumn(0); // depth
param_col[table_num] = enkrvd.getColumn(itab); // itab=[1-7]: krw krg kro krwr krgr krorw krorg
const auto &enkrvdTable = enkrvdTables[table_num];
depth_col[table_num] = enkrvdTable.getDepthColumn();
param_col[table_num] = enkrvdTable.getColumn(itab); // itab=[1-7]: krw krg kro krwr krgr krorw krorg
}
}
} else if (useKeyword && (keyword[0] == 'P' || keyword[1] == 'P') ) {

25
opm/core/simulator/initStateEquil.hpp
View File

@ -27,7 +27,6 @@
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/parser/eclipse/Utility/EquilWrapper.hpp>
#include <opm/parser/eclipse/Utility/SingleRecordTable.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <array>
@ -279,15 +278,15 @@ namespace Opm
// Create Rs functions.
rs_func_.reserve(rec.size());
if (deck->hasKeyword("DISGAS")) {
const std::vector<RsvdTable>& rsvdTables = eclipseState->getRsvdTables();
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].live_oil_table_index > 0) {
if (deck->hasKeyword("RSVD") &&
size_t(rec[i].live_oil_table_index) <= deck->getKeyword("RSVD")->size()) {
Opm::SingleRecordTable rsvd(deck->getKeyword("RSVD"),std::vector<std::string>{"vd", "rs"},rec[i].live_oil_table_index-1);
std::vector<double> vd(rsvd.getColumn("vd"));
std::vector<double> rs(rsvd.getColumn("rs"));
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, vd, rs));
if (rsvdTables.size() > 0 && size_t(rec[i].live_oil_table_index) <= rsvdTables.size()) {
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props,
cell,
rsvdTables[i].getDepthColumn(),
rsvdTables[i].getRsColumn()));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
}
@ -310,15 +309,15 @@ namespace Opm
rv_func_.reserve(rec.size());
if (deck->hasKeyword("VAPOIL")) {
const std::vector<RvvdTable>& rvvdTables = eclipseState->getRvvdTables();
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].wet_gas_table_index > 0) {
if (deck->hasKeyword("RVVD") &&
size_t(rec[i].wet_gas_table_index) <= deck->getKeyword("RVVD")->size()) {
Opm::SingleRecordTable rvvd(deck->getKeyword("RVVD"),std::vector<std::string>{"vd", "rv"},rec[i].wet_gas_table_index-1);
std::vector<double> vd(rvvd.getColumn("vd"));
std::vector<double> rv(rvvd.getColumn("rv"));
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, vd, rv));
if (rvvdTables.size() > 0 && size_t(rec[i].wet_gas_table_index) <= rvvdTables.size()) {
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props,
cell,
rvvdTables[i].getDepthColumn(),
rvvdTables[i].getRvColumn()));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
}

14
tests/liveoil.DATA
View File

@ -1,20 +1,25 @@
RUNSPEC
TABDIMS
-- use the defaults of TABDIMS but the keyword must be present in the deck
-- for it to be usable
/
DIMENS
3 3 3 /
OIL
GAS
WATER
FIELD
GRID
-- tests for the PVT functions need a grid because the OPM-API for the
-- PVT functions assumes the presence of compressed grid cells,
-- i.e. we need to be able to map from compressed to logical cartesian
-- cell indices to find out the PVT region of a cell
DIMENS
3 3 3 /
DXV
1 2 3 /
DYV
@ -22,6 +27,11 @@ DYV
DZV
1 2 3 /
TOPS
9*123.456 /
PROPS
PVTO
-- Rs Pbub Bo Vo
.0 14.7 1.0000 1.20 /

42
tests/test_blackoilfluid.cpp
View File

@ -11,13 +11,9 @@
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
#include <opm/parser/eclipse/Utility/PvtgTable.hpp>
#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
#include <opm/parser/eclipse/Utility/PvcdoTable.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
@ -36,7 +32,7 @@ using namespace Opm;
using namespace std;
std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, PhaseUsage phase_usage_){
std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclipseState, PhaseUsage phase_usage_){
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
int samples = 0;
@ -59,19 +55,20 @@ std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, Pha
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
if (deck->hasKeyword("PVDO")) {
Opm::DeckKeywordConstPtr pvdoKeyword(deck->getKeyword("PVDO"));
const auto& pvdoTables = eclipseState->getPvdoTables();
const auto& pvtoTables = eclipseState->getPvtoTables();
if (!pvdoTables.empty()) {
if (samples > 0) {
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromOil(pvdoKeyword, samples);
splinePvt->initFromOil(pvdoTables, samples);
props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
} else {
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromOil(pvdoKeyword);
deadPvt->initFromOil(pvdoTables);
props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
}
} else if (deck->hasKeyword("PVTO")) {
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO")));
} else if (!pvtoTables.empty()) {
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(pvtoTables));
} else if (deck->hasKeyword("PVCDO")) {
std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
pvcdo->initFromOil(deck->getKeyword("PVCDO"));
@ -83,19 +80,20 @@ std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, Pha
}
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
if (deck->hasKeyword("PVDG")) {
Opm::DeckKeywordConstPtr pvdgKeyword(deck->getKeyword("PVDG"));
const auto& pvdgTables = eclipseState->getPvdgTables();
const auto& pvtgTables = eclipseState->getPvtgTables();
if (!pvdgTables.empty()) {
if (samples > 0) {
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromGas(pvdgKeyword, samples);
splinePvt->initFromGas(pvdgTables, samples);
props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
} else {
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromGas(pvdgKeyword);
deadPvt->initFromGas(pvdgTables);
props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
}
} else if (deck->hasKeyword("PVTG")) {
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG")));
} else if (!pvtgTables.empty()) {
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(pvtgTables));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
}
@ -236,12 +234,11 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
cout << "Reading deck: " << filename << endl;
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr deck(parser->parseFile(filename));
Opm::EclipseStateConstPtr eclipseState(new EclipseState(deck));
// setup pvt interface
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck,phase_usage_);
std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck, eclipseState, phase_usage_);
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
@ -314,10 +311,11 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
cout << "Reading deck: " << filename << endl;
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr deck(parser->parseFile(filename));
Opm::EclipseStateConstPtr eclipseState(new EclipseState(deck));
// setup pvt interface
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck,phase_usage_);
std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck,eclipseState,phase_usage_);
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference

14
tests/wetgas.DATA
View File

@ -1,20 +1,25 @@
RUNSPEC
TABDIMS
-- use the defaults of TABDIMS but the keyword must be present in the deck
-- for it to be usable
/
DIMENS
3 3 3 /
WATER
OIL
GAS
FIELD
GRID
-- tests for the PVT functions need a grid because the OPM-API for the
-- PVT functions assumes the presence of compressed grid cells,
-- i.e. we need to be able to map from compressed to logical cartesian
-- cell indices to find out the PVT region of a cell
DIMENS
3 3 3 /
DXV
1 2 3 /
DYV
@ -22,6 +27,11 @@ DYV
DZV
1 2 3 /
TOPS
9*123.456 /
PROPS
-- PVT PROPERTIES OF DRY GAS (NO VAPOURISED OIL)
-- FROM SPE3 Blackoil Kleppe
--