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

Merge pull request #823 from andlaus/use_fluidstates_for_satfuncs

Browse Source 
Use fluidstates for satfuncs
This commit is contained in:
Atgeirr Flø Rasmussen
2015-07-01 12:57:55 +02:00
parent a127db5638 8ccf1b3dbe
commit 038fb6913d
16 changed files with 846 additions and 476 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
CMakeLists_files.cmake
View File

@@ -98,11 +98,8 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/props/rock/RockBasic.cpp
opm/core/props/rock/RockCompressibility.cpp
opm/core/props/rock/RockFromDeck.cpp
opm/core/props/satfunc/SatFuncGwseg.cpp
opm/core/props/satfunc/SatFuncSimple.cpp
opm/core/props/satfunc/SatFuncStone2.cpp
opm/core/props/satfunc/SaturationPropsBasic.cpp
opm/core/props/satfunc/SaturationPropsFromDeck.cpp
opm/core/simulator/AdaptiveSimulatorTimer.cpp
opm/core/simulator/BlackoilState.cpp
opm/core/simulator/TimeStepControl.cpp
@@ -365,6 +362,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/props/rock/RockFromDeck.hpp
opm/core/props/satfunc/SatFuncBase.hpp
opm/core/props/satfunc/SatFuncGwseg.hpp
opm/core/props/satfunc/SatFuncMultiplexer.hpp
opm/core/props/satfunc/SatFuncSimple.hpp
opm/core/props/satfunc/SatFuncStone2.hpp
opm/core/props/satfunc/SaturationPropsBasic.hpp
@@ -375,7 +373,9 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/simulator/AdaptiveTimeStepping.hpp
opm/core/simulator/AdaptiveTimeStepping_impl.hpp
opm/core/simulator/BlackoilState.hpp
opm/core/simulator/BlackoilStateToFluidState.hpp
opm/core/simulator/EquilibrationHelpers.hpp
opm/core/simulator/ExplicitArraysFluidState.hpp
opm/core/simulator/TimeStepControl.hpp
opm/core/simulator/SimulatorCompressibleTwophase.hpp
opm/core/simulator/SimulatorIncompTwophase.hpp

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

@@ -54,11 +54,10 @@ namespace Opm
rock_.init(eclState, number_of_cells, global_cell, cart_dims);
}
pvt_.init(deck, eclState, /*numSamples=*/0);
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
SaturationPropsFromDeck* ptr
= new SaturationPropsFromDeck();
satprops_.reset(ptr);
ptr->init(deck, eclState, number_of_cells, global_cell, begin_cell_centroids, dimension,
/*numSamples=*/0);
ptr->init(deck, eclState, number_of_cells, global_cell, begin_cell_centroids, dimension);
if (pvt_.numPhases() != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
@@ -89,56 +88,15 @@ namespace Opm
pvt_.init(deck, eclState, pvt_samples);
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", -1);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "gwseg");
if (deck->hasKeyword("ENDSCALE") && threephase_model != "gwseg") {
OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'gwseg' model only.");
}
if (sat_samples > 1) {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
satprops_.reset(ptr);
ptr->init(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegUniform>();
satprops_.reset(ptr);
ptr->init(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else {
OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
}
} else {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
satprops_.reset(ptr);
ptr->init(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else {
OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
}
}
SaturationPropsFromDeck* ptr
= new SaturationPropsFromDeck();
satprops_.reset(ptr);
ptr->init(deck, eclState, number_of_cells, global_cell, begin_cell_centroids, dimension);
if (pvt_.numPhases() != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("

2
opm/core/props/IncompPropertiesFromDeck.cpp
View File

@@ -32,7 +32,7 @@ namespace Opm
{
rock_.init(eclState, grid.number_of_cells, grid.global_cell, grid.cartdims);
pvt_.init(deck);
satprops_.init(deck, eclState, grid, 200);
satprops_.init(deck, eclState, grid);
if (pvt_.numPhases() != satprops_.numPhases()) {
OPM_THROW(std::runtime_error, "IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");

2
opm/core/props/IncompPropertiesFromDeck.hpp
View File

@@ -138,7 +138,7 @@ namespace Opm
private:
RockFromDeck rock_;
PvtPropertiesIncompFromDeck pvt_;
SaturationPropsFromDeck<SatFuncStone2Uniform> satprops_;
SaturationPropsFromDeck satprops_;
};

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

@@ -77,6 +77,8 @@ namespace Opm
class SatFuncBase : public BlackoilPhases
{
public:
virtual ~SatFuncBase() {};
void init(Opm::EclipseStateConstPtr eclipseState,
const int table_num,
const PhaseUsage phase_usg,

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

@@ -1,33 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/props/satfunc/SatFuncGwseg.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>
namespace Opm
{
// SatFuncGwseg.cpp can now be removed
} // namespace Opm

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

@@ -27,34 +27,42 @@ namespace Opm
class SatFuncGwseg : public SatFuncBase<TableType>
{
public:
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 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:
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr) const;
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const;
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc) const;
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const;
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst) const;
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const;
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const;
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const;
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc, const EPSTransforms* epst) const;
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds, const EPSTransforms* epst) const;
};
typedef SatFuncGwseg<UniformTableLinear<double> > SatFuncGwsegUniform;
typedef SatFuncGwseg<NonuniformTableLinear<double> > SatFuncGwsegNonuniform;
template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, double* kr) const
template <class FluidState>
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr) const
{
if (this->phase_usage.num_phases == 3) {
// 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 sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-5;
swco = std::min(swco,sw-eps);
@@ -77,7 +85,7 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
@@ -86,7 +94,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
@@ -96,7 +104,8 @@ namespace Opm
template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, double* kr, const EPSTransforms* epst) const
template <class FluidState>
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst) const
{
if (this->phase_usage.num_phases == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
@@ -108,8 +117,8 @@ namespace Opm
// 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 sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
@@ -141,7 +150,7 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
@@ -150,7 +159,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
@@ -159,7 +168,8 @@ namespace Opm
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
template <class FluidState>
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{
if (this->phase_usage.num_phases == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
@@ -171,8 +181,8 @@ namespace Opm
// 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 sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
@@ -233,7 +243,7 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
@@ -242,7 +252,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
@@ -251,7 +261,8 @@ namespace Opm
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds) const
template<class FluidState>
void SatFuncGwseg<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
@@ -261,8 +272,8 @@ namespace Opm
// and gas, with oil present in both water and gas zones.
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 sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-5;
swco = std::min(swco,sw-eps);
@@ -297,7 +308,7 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
@@ -310,7 +321,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
@@ -324,7 +335,8 @@ namespace Opm
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const
template<class FluidState>
void SatFuncGwseg<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
@@ -339,8 +351,8 @@ namespace Opm
// 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 sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
@@ -386,7 +398,7 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
@@ -399,7 +411,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
@@ -412,7 +424,8 @@ namespace Opm
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
template<class FluidState>
void SatFuncGwseg<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
@@ -427,8 +440,8 @@ namespace Opm
// 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 sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
@@ -506,7 +519,7 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
@@ -519,7 +532,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
@@ -532,38 +545,41 @@ namespace Opm
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPc(const double* s, double* pc) const
template <class FluidState>
void SatFuncGwseg<TableType>::evalPc(const FluidState& fluidState, double* pc) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]);
pc[pos] = this->pcow_(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
pc[pos] = this->pcog_(fluidState.saturation(pos));
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPc(const double* s, double* pc, const EPSTransforms* epst) const
template <class FluidState>
void SatFuncGwseg<TableType>::evalPc(const FluidState& fluidState, double* pc, const EPSTransforms* epst) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
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]);
double _sw = epst->wat.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->wat.pcFactor*this->pcow_(_sw);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
double _sg = epst->gas.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]);
double _sg = epst->gas.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->gas.pcFactor*this->pcog_(_sg);
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPcDeriv(const double* s, double* pc, double* dpcds) const
template <class FluidState>
void SatFuncGwseg<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
@@ -576,19 +592,20 @@ namespace Opm
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
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]);
pc[pos] = this->pcow_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcow_.derivative(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]);
pc[pos] = this->pcog_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcog_.derivative(fluidState.saturation(pos));
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const
template <class FluidState>
void SatFuncGwseg<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds, const EPSTransforms* epst) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
@@ -601,16 +618,16 @@ namespace Opm
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
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]);
double _sw = epst->wat.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->wat.pcFactor*this->pcow_(_sw);
double _dsdsw = epst->wat.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]);
double _dsdsw = epst->wat.scaleSatDerivPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
dpcds[np*pos + pos] = epst->wat.pcFactor*_dsdsw*this->pcow_.derivative(_sw);
}
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]);
double _sg = epst->gas.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
pc[pos] = epst->gas.pcFactor*this->pcog_(_sg);
double _dsdsg = epst->gas.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]);
double _dsdsg = epst->gas.scaleSatDerivPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
dpcds[np*pos + pos] = epst->gas.pcFactor*_dsdsg*this->pcog_.derivative(_sg);
}
}

261
opm/core/props/satfunc/SatFuncMultiplexer.hpp Normal file
View File

@@ -0,0 +1,261 @@
/*
Copyright 2015 Andreas Lauser
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
/*!
* \file
* \copydoc Opm::SatFuncMultiplexer
*/
#ifndef OPM_SAT_FUNC_MULTIPLEXER_HPP
#define OPM_SAT_FUNC_MULTIPLEXER_HPP
#include "SatFuncGwseg.hpp"
#include "SatFuncSimple.hpp"
#include "SatFuncStone2.hpp"
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <memory>
namespace Opm
{
/*!
* \brief This is a multiplexer class which allows to saturation function to be used at
* runtime.
*/
class SatFuncMultiplexer
{
public:
typedef NonuniformTableLinear<double> TableType;
enum SatFuncType {
Gwseg, // <- supposed to be the default
Stone2,
Simple
};
// this is a helper macro which helps to save us from RSI in the following
#define OPM_MULTIPLEXER_CONST const
#define OPM_MULTIPLEXER_NON_CONST
#define OPM_MULTIPLEXER_SATFUNC_CALL(codeToCall, CONST) \
switch (satFuncType_) { \
case Gwseg: { \
typedef SatFuncGwseg<TableType> SatFunc; \
__attribute__((unused)) CONST SatFunc& satFunc = \
*static_cast<SatFunc*>(satFunc_.get()); \
codeToCall; \
break; \
} \
case Stone2: { \
typedef SatFuncStone2<TableType> SatFunc; \
__attribute__((unused)) CONST SatFunc& satFunc = \
*static_cast<SatFunc*>(satFunc_.get()); \
codeToCall; \
break; \
} \
case Simple: { \
typedef SatFuncSimple<TableType> SatFunc; \
__attribute__((unused)) CONST SatFunc& satFunc = \
*static_cast<SatFunc*>(satFunc_.get()); \
codeToCall; \
break; \
} \
};
SatFuncMultiplexer()
{}
SatFuncMultiplexer(const SatFuncMultiplexer&)
{}
~SatFuncMultiplexer()
{}
/*!
* \brief Pick the correct saturation function type and initialize the object using
* an ECL deck.
*/
void initFromDeck(Opm::EclipseStateConstPtr eclState,
size_t tableIdx,
SatFuncType satFuncType)
{
auto phaseUsage = phaseUsageFromDeck(eclState);
setSatFuncType(satFuncType);
OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.init(eclState,
tableIdx,
phaseUsage,
/*samples=*/0), OPM_MULTIPLEXER_NON_CONST);
}
/*!
* \brief Returns the type of the actually selected saturation function.
*/
void setSatFuncType(SatFuncType newSatFuncType)
{
// set the type of the saturation function
satFuncType_ = newSatFuncType;
// call the default constructor for the selected saturation function
switch (satFuncType_) {
case Gwseg: {
typedef SatFuncGwseg<TableType> SatFunc;
satFunc_.reset(new SatFunc);
break;
}
case Stone2: {
typedef SatFuncStone2<TableType> SatFunc;
satFunc_.reset(new SatFunc);
break;
}
case Simple: {
typedef SatFuncSimple<TableType> SatFunc;
satFunc_.reset(new SatFunc);
break;
}
}
}
/*!
* \brief Returns which saturation function is actually selected.
*/
SatFuncType satFuncType() const
{ return satFuncType_; }
/*!
* \brief Returns a pointer to the base saturation function.
*
* All underlying saturation functions derive from this class. It stores the stuff
* which is common for all saturation functions (e.g. the endpoint scaling
* parameters).
*/
const SatFuncBase<TableType>& getSatFuncBase() const
{ return *satFunc_; }
SatFuncBase<TableType>& getSatFuncBase()
{ return *satFunc_; }
/*!
* \brief Return the raw "Gwseg" saturation function object.
*
* If the saturation function type is not "Gwseg", an assertation is triggered for
* debug builds.
*/
SatFuncGwseg<TableType>& getGwseg()
{
assert(satFuncType_ == Gwseg);
return *static_cast<SatFuncGwseg<TableType>*>(satFunc_.get());
}
const SatFuncGwseg<TableType>& getGwseg() const
{
assert(satFuncType_ == Gwseg);
return *static_cast<const SatFuncGwseg<TableType>*>(satFunc_.get());
}
/*!
* \brief Return the raw "Simple" saturation function object.
*
* If the saturation function type is not "Simple", an assertation is triggered for
* debug builds.
*/
SatFuncSimple<TableType>& getSimple()
{
assert(satFuncType_ == Simple);
return *static_cast<SatFuncSimple<TableType>*>(satFunc_.get());
}
const SatFuncSimple<TableType>& getSimple() const
{
assert(satFuncType_ == Simple);
return *static_cast<const SatFuncSimple<TableType>*>(satFunc_.get());
}
/*!
* \brief Return the raw "Stone2" saturation function object.
*
* If the saturation function type is not "Stone2", an assertation is triggered for
* debug builds.
*/
SatFuncStone2<TableType>& getStone2()
{
assert(satFuncType_ == Stone2);
return *static_cast<SatFuncStone2<TableType>*>(satFunc_.get());
}
const SatFuncStone2<TableType>& getStone2() const
{
assert(satFuncType_ == Stone2);
return *static_cast<const SatFuncStone2<TableType>*>(satFunc_.get());
}
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKr(fluidState, kr), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKrDeriv(fluidState, kr, dkrds), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPc(fluidState, pc), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPcDeriv(fluidState, pc, dpcds), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKr(fluidState, kr, epst), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKr(fluidState, kr, epst, epst_hyst, sat_hyst), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKrDeriv(fluidState, kr, dkrds, epst), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKrDeriv(fluidState, kr, dkrds, epst, epst_hyst, sat_hyst), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc, const EPSTransforms* epst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPc(fluidState, pc, epst), OPM_MULTIPLEXER_CONST); }
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds, const EPSTransforms* epst) const
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPcDeriv(fluidState, pc, dpcds, epst), OPM_MULTIPLEXER_CONST); }
// undefine the helper macro here because we don't need it anymore
#undef OPM_MULTIPLEXER_SATFUNC_CALL
#undef OPM_MULTIPLEXER_NON_CONST
#undef OPM_MULTIPLEXER_CONST
private:
SatFuncType satFuncType_;
std::unique_ptr<SatFuncBase<TableType> > satFunc_;
};
} // namespace Opm
#endif // OPM_SAT_FUNC_MULTIPLEXER_HPP

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

@@ -27,21 +27,31 @@ namespace Opm
class SatFuncSimple : public SatFuncBase<TableType>
{
public:
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;
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr) const;
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const;
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc) const;
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const;
void evalKr(const double* /* s */, double* /* kr */, const EPSTransforms* /* epst */) const
template <class FluidState>
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
void evalKr(const double* /* s */, double* /* kr */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
template <class FluidState>
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
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
template <class FluidState>
void evalKrDeriv(const FluidState& /* fluidState */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */) const;
template <class FluidState>
void evalKrDeriv(const FluidState& /* fluidState */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
void evalPc(const double* /* s */, double* /* pc */, const EPSTransforms* /* epst */) const
template <class FluidState>
void evalPc(const FluidState& /* fluidState */, 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
template <class FluidState>
void evalPcDeriv(const FluidState& /* fluidState */, double* /* pc */, double* /* dpcds */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
private:
@@ -52,12 +62,13 @@ namespace Opm
typedef SatFuncSimple<NonuniformTableLinear<double> > SatFuncSimpleNonuniform;
template<class TableType>
void SatFuncSimple<TableType>::evalKr(const double* s, double* kr) const
template<class FluidState>
void SatFuncSimple<TableType>::evalKr(const FluidState& fluidState, double* kr) const
{
if (this->phase_usage.num_phases == 3) {
// A simplified relative permeability model.
double sw = s[BlackoilPhases::Aqua];
double sg = s[BlackoilPhases::Vapour];
double sw = fluidState.saturation(BlackoilPhases::Aqua);
double sg = fluidState.saturation(BlackoilPhases::Vapour);
double krw = this->krw_(sw);
double krg = this->krg_(sg);
double krow = this->krow_(sw + sg); // = 1 - so
@@ -75,9 +86,9 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double so = s[opos];
double so = fluidState.saturation(opos);
double krow = this->krow_(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
@@ -85,7 +96,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
@@ -94,15 +105,16 @@ namespace Opm
}
template<class TableType>
void SatFuncSimple<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds) const
template<class FluidState>
void SatFuncSimple<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// A simplified relative permeability model.
double sw = s[BlackoilPhases::Aqua];
double sg = s[BlackoilPhases::Vapour];
double sw = fluidState.saturation(BlackoilPhases::Aqua);
double sg = fluidState.saturation(BlackoilPhases::Vapour);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krg = this->krg_(sg);
@@ -133,10 +145,10 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double so = s[opos];
double so = fluidState.saturation(opos);
double krow = this->krow_(1.0-so);
double dkrow = this->krow_.derivative(1.0-so);
kr[wpos] = krw;
@@ -147,7 +159,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
@@ -161,7 +173,8 @@ namespace Opm
}
template<class TableType>
void SatFuncSimple<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const
template<class FluidState>
void SatFuncSimple<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
@@ -172,24 +185,24 @@ namespace Opm
// 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 _sw = epst->wat.scaleSat(fluidState.saturation(wpos), 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _dsdsw = epst->wat.scaleSatDeriv(fluidState.saturation(wpos), 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _sg = epst->gas.scaleSat(fluidState.saturation(gpos), 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _dsdsg = epst->gas.scaleSatDeriv(fluidState.saturation(gpos), 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _sow = epst->watoil.scaleSat(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdsow = epst->watoil.scaleSatDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
//double _sog = epst->gasoil.scaleSat(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
//double _dsdsog = epst->gasoil.scaleSatDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double krw = epst->wat.scaleKr(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));
double krw = epst->wat.scaleKr(fluidState.saturation(wpos), this->krw_(_sw), this->krwr_);
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(fluidState.saturation(wpos), this->krw_.derivative(_sw));
double krg = epst->gas.scaleKr(fluidState.saturation(gpos), this->krg_(_sg), this->krgr_);
double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(fluidState.saturation(gpos), this->krg_.derivative(_sg));
// TODO Check the arguments to the krow- and krog-tables below...
double krow = epst->watoil.scaleKr(1.0-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 krow = epst->watoil.scaleKr(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krow_(1.0-_sow-this->smin_[gpos]), this->krorw_); // ????
double dkrow = _dsdsow*epst->watoil.scaleKrDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krow_.derivative(1.0-_sow-this->smin_[gpos])); // ????
//double krog = epst->gasoil.scaleKr(this->krog_(1.0-_sog-this->smin_[wpos]), 1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krorg_); // ????
//double dkrog = _dsdsog*epst->gasoil.scaleKrDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krog_.derivative(1.0-_sog-this->smin_[wpos])); // ????
// double krocw = krocw_;
kr[wpos] = krw;
kr[gpos] = krg;
@@ -213,10 +226,10 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double so = s[opos];
double so = fluidState.saturation(opos);
double krow = this->krow_(1.0-so);
double dkrow = this->krow_.derivative(1.0-so);
kr[wpos] = krw;
@@ -227,7 +240,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
@@ -241,21 +254,23 @@ namespace Opm
}
template<class TableType>
void SatFuncSimple<TableType>::evalPc(const double* s, double* pc) const
template<class FluidState>
void SatFuncSimple<TableType>::evalPc(const FluidState& fluidState, double* pc) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]);
pc[pos] = this->pcow_(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
pc[pos] = this->pcog_(fluidState.saturation(pos));
}
}
template<class TableType>
void SatFuncSimple<TableType>::evalPcDeriv(const double* s, double* pc, double* dpcds) const
template<class FluidState>
void SatFuncSimple<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
@@ -268,13 +283,13 @@ namespace Opm
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]);
pc[pos] = this->pcow_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcow_.derivative(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]);
pc[pos] = this->pcog_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcog_.derivative(fluidState.saturation(pos));
}
}

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

@@ -1,33 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/props/satfunc/SatFuncStone2.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>
namespace Opm
{
// SatFuncStone2.cpp can now be removed
} // namespace Opm

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

@@ -27,22 +27,32 @@ namespace Opm
class SatFuncStone2 : public SatFuncBase<TableType>
{
public:
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;
template <class FluidState>
void evalKr(const FluidState& fluidState, double* kr) const;
template <class FluidState>
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const;
template <class FluidState>
void evalPc(const FluidState& fluidState, double* pc) const;
template <class FluidState>
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const;
void evalKr(const double* /* s */, double* /* kr */, const EPSTransforms* /* epst */) const
template <class FluidState>
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalKr(const double* /* s */, double* /* kr */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
template <class FluidState>
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalKrDeriv(const double* /* s */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */) const
template <class FluidState>
void evalKrDeriv(const FluidState& /* fluidState */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalKrDeriv(const double* /* s */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
template <class FluidState>
void evalKrDeriv(const FluidState& /* fluidState */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */, const 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
template <class FluidState>
void evalPc(const FluidState& /* fluidState */, 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
template <class FluidState>
void evalPcDeriv(const FluidState& /* fluidState */, double* /* pc */, double* /* dpcds */, const EPSTransforms* /* epst */) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
private:
@@ -53,12 +63,13 @@ namespace Opm
typedef SatFuncStone2<NonuniformTableLinear<double> > SatFuncStone2Nonuniform;
template<class TableType>
void SatFuncStone2<TableType>::evalKr(const double* s, double* kr) const
template <class FluidState>
void SatFuncStone2<TableType>::evalKr(const FluidState& fluidState, double* kr) const
{
if (this->phase_usage.num_phases == 3) {
// Stone-II relative permeability model.
double sw = s[BlackoilPhases::Aqua];
double sg = s[BlackoilPhases::Vapour];
double sw = fluidState.saturation(BlackoilPhases::Aqua);
double sg = fluidState.saturation(BlackoilPhases::Vapour);
double krw = this->krw_(sw);
double krg = this->krg_(sg);
double krow = this->krow_(sw + sg); // = 1 - so
@@ -76,7 +87,7 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
@@ -85,7 +96,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
@@ -94,15 +105,16 @@ namespace Opm
}
template<class TableType>
void SatFuncStone2<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds) const
template <class FluidState>
void SatFuncStone2<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Stone-II relative permeability model.
double sw = s[BlackoilPhases::Aqua];
double sg = s[BlackoilPhases::Vapour];
double sw = fluidState.saturation(BlackoilPhases::Aqua);
double sg = fluidState.saturation(BlackoilPhases::Vapour);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krg = this->krg_(sg);
@@ -129,7 +141,7 @@ namespace Opm
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 sw = fluidState.saturation(wpos);
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
@@ -142,7 +154,7 @@ namespace Opm
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 sg = fluidState.saturation(gpos);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
@@ -156,21 +168,23 @@ namespace Opm
}
template<class TableType>
void SatFuncStone2<TableType>::evalPc(const double* s, double* pc) const
template <class FluidState>
void SatFuncStone2<TableType>::evalPc(const FluidState& fluidState, double* pc) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]);
pc[pos] = this->pcow_(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
pc[pos] = this->pcog_(fluidState.saturation(pos));
}
}
template<class TableType>
void SatFuncStone2<TableType>::evalPcDeriv(const double* s, double* pc, double* dpcds) const
template <class FluidState>
void SatFuncStone2<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
@@ -183,13 +197,13 @@ namespace Opm
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]);
pc[pos] = this->pcow_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcow_.derivative(fluidState.saturation(pos));
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]);
pc[pos] = this->pcog_(fluidState.saturation(pos));
dpcds[np*pos + pos] = this->pcog_.derivative(fluidState.saturation(pos));
}
}

35
opm/core/props/satfunc/SaturationPropsFromDeck.cpp
View File

@@ -1,35 +0,0 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>
namespace Opm
{
// This file should be removed in the future.
// Holding off until refactoring of SaturationPropsFromDeck class is done.
} // namespace Opm

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

@@ -23,9 +23,7 @@
#include <opm/core/props/satfunc/SaturationPropsInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SatFuncStone2.hpp>
#include <opm/core/props/satfunc/SatFuncSimple.hpp>
#include <opm/core/props/satfunc/SatFuncGwseg.hpp>
#include <opm/core/props/satfunc/SatFuncMultiplexer.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@@ -36,33 +34,21 @@ struct UnstructuredGrid;
namespace Opm
{
/// Interface to saturation functions from deck.
/// Possible values for template argument (for now):
/// SatFuncSetStone2Nonuniform,
/// SatFuncSetStone2Uniform.
/// SatFuncSetSimpleNonuniform,
/// SatFuncSetSimpleUniform.
template <class SatFuncSet>
class SaturationPropsFromDeck : public SaturationPropsInterface
{
public:
/// Default constructor.
SaturationPropsFromDeck();
inline SaturationPropsFromDeck();
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
/// \param[in] samples Number of uniform sample points for saturation tables.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
void init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
const UnstructuredGrid& grid,
const int samples);
inline void init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
const UnstructuredGrid& grid);
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
@@ -75,19 +61,16 @@ namespace Opm
/// global cell indices used in the deck.
/// \param[in] begin_cell_centroids Pointer to the first cell_centroid of the grid.
/// \param[in] dimensions The dimensions of the grid.
/// \param[in] samples Number of uniform sample points for saturation tables.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
template<class T>
void init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const int samples);
inline void init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
/// \return P, the number of phases.
int numPhases() const;
inline int numPhases() const;
/// Relative permeability.
/// \param[in] n Number of data points.
@@ -98,11 +81,11 @@ namespace Opm
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const;
inline void relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const;
/// Capillary pressure.
/// \param[in] n Number of data points.
@@ -113,39 +96,48 @@ namespace Opm
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const;
inline void capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const;
/// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void satRange(const int n,
const int* cells,
double* smin,
double* smax) const;
inline void satRange(const int n,
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);
inline void updateSatHyst(const int n,
const int* cells,
const double* s);
/// Update capillary pressure scaling according to pressure diff. and initial water saturation.
/// \param[in] cell Cell index.
/// \param[in] pcow P_oil - P_water.
/// \param[in/out] swat Water saturation. / Possibly modified Water saturation.
void swatInitScaling(const int cell,
const double pcow,
double & swat);
inline void swatInitScaling(const int cell,
const double pcow,
double & swat);
private:
// internal helper method for satRange()
template <class SaturationFunction>
void satRange_(const SaturationFunction& satFunc,
const int cellIdx,
const int* cells,
double* smin,
double* smax) const;
PhaseUsage phase_usage_;
std::vector<SatFuncSet> satfuncset_;
typedef Opm::SatFuncMultiplexer SatFunc;
std::vector<SatFunc> satfunc_;
std::vector<int> cell_to_func_; // = SATNUM - 1
std::vector<int> cell_to_func_imb_;
@@ -156,9 +148,6 @@ namespace Opm
std::vector<EPSTransforms> eps_transf_hyst_;
std::vector<SatHyst> sat_hyst_;
typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
template<class T>
void initEPS(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,

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

@@ -24,6 +24,7 @@
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/simulator/ExplicitArraysFluidState.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridHelpers.hpp>
@@ -41,39 +42,36 @@ namespace Opm
/// Default constructor.
template <class SatFuncSet>
SaturationPropsFromDeck<SatFuncSet>::SaturationPropsFromDeck()
inline
SaturationPropsFromDeck::SaturationPropsFromDeck()
{
}
/// Initialize from deck.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr deck,
inline
void SaturationPropsFromDeck::init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
const UnstructuredGrid& grid,
const int samples)
const UnstructuredGrid& grid)
{
this->init(deck, eclipseState, grid.number_of_cells,
grid.global_cell, grid.cell_centroids,
grid.dimensions, samples);
grid.dimensions);
}
/// Initialize from deck.
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr deck,
void SaturationPropsFromDeck::init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const int samples)
int dimensions)
{
phase_usage_ = phaseUsageFromDeck(deck);
// Extract input data.
// Oil phase should be active.
if (!phase_usage_.phase_used[Liquid]) {
if (!phase_usage_.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- oil phase must be active.");
}
@@ -93,11 +91,11 @@ namespace Opm
// Obtain SATNUM, if it exists, and create cell_to_func_.
// Otherwise, let the cell_to_func_ mapping be just empty.
int satfuncs_expected = 1;
cell_to_func_.resize(number_of_cells, /*value=*/0);
if (deck->hasKeyword("SATNUM")) {
const std::vector<int>& satnum = deck->getKeyword("SATNUM")->getIntData();
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
const int num_cells = number_of_cells;
cell_to_func_.resize(num_cells);
const int* gc = global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int deck_pos = (gc == NULL) ? cell : gc[cell];
@@ -108,13 +106,13 @@ namespace Opm
// Find number of tables, check for consistency.
enum { Uninitialized = -1 };
int num_tables = Uninitialized;
if (phase_usage_.phase_used[Aqua]) {
if (phase_usage_.phase_used[BlackoilPhases::Aqua]) {
num_tables = deck->getKeyword("SWOF")->size();
if (num_tables < satfuncs_expected) {
OPM_THROW(std::runtime_error, "Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
}
}
if (phase_usage_.phase_used[Vapour]) {
if (phase_usage_.phase_used[BlackoilPhases::Vapour]) {
int num_sgof_tables = deck->getKeyword("SGOF")->size();
if (num_sgof_tables < satfuncs_expected) {
OPM_THROW(std::runtime_error, "Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
@@ -126,12 +124,17 @@ namespace Opm
}
}
// Initialize tables.
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(eclipseState, table, phase_usage_, samples);
// Initialize saturation function objects.
satfunc_.resize(num_tables);
SatFuncMultiplexer::SatFuncType satFuncType = SatFuncMultiplexer::Gwseg;
if (deck->hasKeyword("STONE2")) {
satFuncType = SatFuncMultiplexer::Stone2;
}
for (int table = 0; table < num_tables; ++table) {
satfunc_[table].initFromDeck(eclipseState, table, satFuncType);
}
// Check EHYSTR status
do_hyst_ = false;
if (hysteresis_switch && deck->hasKeyword("EHYSTR")) {
@@ -249,8 +252,8 @@ namespace Opm
/// \return P, the number of phases.
template <class SatFuncSet>
int SaturationPropsFromDeck<SatFuncSet>::numPhases() const
inline
int SaturationPropsFromDeck::numPhases() const
{
return phase_usage_.num_phases;
}
@@ -268,8 +271,8 @@ namespace Opm
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dkr_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::relperm(const int n,
inline
void SaturationPropsFromDeck::relperm(const int n,
const double* s,
const int* cells,
double* kr,
@@ -277,27 +280,31 @@ namespace Opm
{
assert(cells != 0);
ExplicitArraysFluidState fluidState;
fluidState.setSaturationArray(s);
const int np = phase_usage_.num_phases;
if (dkrds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
fluidState.setIndex(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]]));
satfunc_[cell_to_func_[cells[i]]].evalKrDeriv(fluidState, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]]));
} else if (do_eps_) {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]]));
satfunc_[cell_to_func_[cells[i]]].evalKrDeriv(fluidState, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]]));
} else {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
satfunc_[cell_to_func_[cells[i]]].evalKrDeriv(fluidState, kr + np*i, dkrds + np*np*i);
}
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++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]]));
satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]]));
} else if (do_eps_) {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i, &(eps_transf_[cells[i]]));
satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i, &(eps_transf_[cells[i]]));
} else {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i);
satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i);
}
}
}
@@ -316,8 +323,8 @@ namespace Opm
/// The P^2 derivative matrix is
/// m_{ij} = \frac{dpc_i}{ds^j},
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::capPress(const int n,
inline
void SaturationPropsFromDeck::capPress(const int n,
const double* s,
const int* cells,
double* pc,
@@ -325,41 +332,67 @@ namespace Opm
{
assert(cells != 0);
ExplicitArraysFluidState fluidState;
fluidState.setSaturationArray(s);
const int np = phase_usage_.num_phases;
if (dpcds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
fluidState.setIndex(i);
if (do_eps_) {
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i, &(eps_transf_[cells[i]]));
satfunc_[cell_to_func_[cells[i]]].evalPcDeriv(fluidState, pc + np*i, dpcds + np*np*i, &(eps_transf_[cells[i]]));
} else {
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i);
satfunc_[cell_to_func_[cells[i]]].evalPcDeriv(fluidState, pc + np*i, dpcds + np*np*i);
}
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
fluidState.setIndex(i);
if (do_eps_) {
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i, &(eps_transf_[cells[i]]));
satfunc_[cell_to_func_[cells[i]]].evalPc(fluidState, pc + np*i, &(eps_transf_[cells[i]]));
} else {
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i);
satfunc_[cell_to_func_[cells[i]]].evalPc(fluidState, pc + np*i);
}
}
}
}
/// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::satRange(const int n,
inline
void SaturationPropsFromDeck::satRange(const int n,
const int* cells,
double* smin,
double* smax) const
{
for (int cellIdx = 0; cellIdx < n; ++cellIdx) {
const SatFuncMultiplexer& multiplexerSatFunc = satfunc_[cell_to_func_[cellIdx]];
switch (multiplexerSatFunc.satFuncType()) {
case SatFuncMultiplexer::Gwseg:
satRange_(multiplexerSatFunc.getGwseg(), cellIdx, cells, smin, smax);
break;
case SatFuncMultiplexer::Stone2:
satRange_(multiplexerSatFunc.getStone2(), cellIdx, cells, smin, smax);
break;
case SatFuncMultiplexer::Simple:
satRange_(multiplexerSatFunc.getSimple(), cellIdx, cells, smin, smax);
break;
}
}
}
template <class SaturationFunction>
void SaturationPropsFromDeck::satRange_(const SaturationFunction& satFunc,
const int cellIdx,
const int* cells,
double* smin,
double* smax) const
{
assert(cells != 0);
const int np = phase_usage_.num_phases;
@@ -368,35 +401,32 @@ namespace Opm
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 ? funcForCell(cells[i]).smin_[wpos]
: eps_transf_[cells[i]].wat.smin;
smax[np*i + wpos] = eps_transf_[cells[i]].wat.doNotScale ? funcForCell(cells[i]).smax_[wpos]
: eps_transf_[cells[i]].wat.smax;
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]].gas.doNotScale ? funcForCell(cells[i]).smin_[gpos]
: eps_transf_[cells[i]].gas.smin;
smax[np*i + gpos] = eps_transf_[cells[i]].gas.doNotScale ? funcForCell(cells[i]).smax_[gpos]
: eps_transf_[cells[i]].gas.smax;
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]);
}
smin[np*cellIdx + opos] = 1.0;
smax[np*cellIdx + opos] = 1.0;
if (phase_usage_.phase_used[BlackoilPhases::Aqua]) {
smin[np*cellIdx + wpos] = eps_transf_[cells[cellIdx]].wat.doNotScale ? satFunc.smin_[wpos]
: eps_transf_[cells[cellIdx]].wat.smin;
smax[np*cellIdx + wpos] = eps_transf_[cells[cellIdx]].wat.doNotScale ? satFunc.smax_[wpos]
: eps_transf_[cells[cellIdx]].wat.smax;
smin[np*cellIdx + opos] -= smax[np*cellIdx + wpos];
smax[np*cellIdx + opos] -= smin[np*cellIdx + wpos];
}
if (phase_usage_.phase_used[BlackoilPhases::Vapour]) {
smin[np*cellIdx + gpos] = eps_transf_[cells[cellIdx]].gas.doNotScale ? satFunc.smin_[gpos]
: eps_transf_[cells[cellIdx]].gas.smin;
smax[np*cellIdx + gpos] = eps_transf_[cells[cellIdx]].gas.doNotScale ? satFunc.smax_[gpos]
: eps_transf_[cells[cellIdx]].gas.smax;
smin[np*cellIdx + opos] -= smax[np*cellIdx + gpos];
smax[np*cellIdx + opos] -= smin[np*cellIdx + gpos];
}
if (phase_usage_.phase_used[BlackoilPhases::Vapour] && phase_usage_.phase_used[BlackoilPhases::Aqua]) {
smin[np*cellIdx + opos] = std::max(0.0,smin[np*cellIdx + 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];
}
for (int p = 0; p < np; ++p) {
smin[np*cellIdx + p] = satFunc.smin_[p];
smax[np*cellIdx + p] = satFunc.smax_[p];
}
}
}
@@ -405,8 +435,8 @@ namespace Opm
/// 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,
inline
void SaturationPropsFromDeck::updateSatHyst(const int n,
const int* cells,
const double* s)
{
@@ -416,7 +446,7 @@ namespace Opm
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]]));
satfunc_[cell_to_func_[cells[i]]].getSatFuncBase().updateSatHyst(s + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]]));
}
}
}
@@ -426,8 +456,8 @@ namespace Opm
/// \param[in] cell Cell index.
/// \param[in] pcow P_oil - P_water.
/// \param[in/out] swat Water saturation. / Possibly modified Water saturation.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::swatInitScaling(const int cell,
inline
void SaturationPropsFromDeck::swatInitScaling(const int cell,
const double pcow,
double& swat)
{
@@ -443,8 +473,11 @@ namespace Opm
const int max_np = BlackoilPhases::MaxNumPhases;
double s[max_np] = { 0.0 };
s[wpos] = swat;
ExplicitArraysFluidState fluidState;
fluidState.setSaturationArray(s);
fluidState.setIndex(0);
double pc[max_np] = { 0.0 };
funcForCell(cell).evalPc(s, pc, &(eps_transf_[cell]));
satfunc_[cell_to_func_[cell]].evalPc(fluidState, pc, &(eps_transf_[cell]));
if (pc[wpos] > pc_low_threshold) {
eps_transf_[cell].wat.pcFactor *= pcow/pc[wpos];
}
@@ -455,18 +488,9 @@ namespace Opm
}
// Map the cell number to the correct function set.
template <class SatFuncSet>
const typename SaturationPropsFromDeck<SatFuncSet>::Funcs&
SaturationPropsFromDeck<SatFuncSet>::funcForCell(const int cell) const
{
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
// Initialize saturation scaling parameters
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(Opm::DeckConstPtr deck,
void SaturationPropsFromDeck::initEPS(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
@@ -485,56 +509,57 @@ namespace Opm
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];
const bool oilWater = phase_usage_.phase_used[BlackoilPhases::Aqua] && phase_usage_.phase_used[BlackoilPhases::Liquid] && !phase_usage_.phase_used[BlackoilPhases::Vapour];
const bool oilGas = !phase_usage_.phase_used[BlackoilPhases::Aqua] && phase_usage_.phase_used[BlackoilPhases::Liquid] && phase_usage_.phase_used[BlackoilPhases::Vapour];
const bool threephase = phase_usage_.phase_used[BlackoilPhases::Aqua] && phase_usage_.phase_used[BlackoilPhases::Liquid] && phase_usage_.phase_used[BlackoilPhases::Vapour];
for (int cell = 0; cell < number_of_cells; ++cell) {
auto& satFuncBase = satfunc_[cell_to_func_[cell]].getSatFuncBase();
if (threephase || oilWater) {
// ### krw
initEPSParam(cell, eps_transf[cell].wat, false,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
funcForCell(cell).smax_[wpos],
funcForCell(cell).sowcr_,
oilWater ? -1.0 : funcForCell(cell).smin_[gpos],
funcForCell(cell).krwr_,
funcForCell(cell).krwmax_,
funcForCell(cell).pcwmax_,
satFuncBase.smin_[wpos],
satFuncBase.swcr_,
satFuncBase.smax_[wpos],
satFuncBase.sowcr_,
oilWater ? -1.0 : satFuncBase.smin_[gpos],
satFuncBase.krwr_,
satFuncBase.krwmax_,
satFuncBase.pcwmax_,
eps_vec[0], eps_vec[2], eps_vec[1], eps_vec[6], eps_vec[3], eps_vec[11], eps_vec[8], eps_vec[15]);
// ### krow
initEPSParam(cell, eps_transf[cell].watoil, true,
0.0,
funcForCell(cell).sowcr_,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
oilWater ? -1.0 : funcForCell(cell).smin_[gpos],
funcForCell(cell).krorw_,
funcForCell(cell).kromax_,
satFuncBase.sowcr_,
satFuncBase.smin_[wpos],
satFuncBase.swcr_,
oilWater ? -1.0 : satFuncBase.smin_[gpos],
satFuncBase.krorw_,
satFuncBase.kromax_,
0.0,
eps_vec[0], eps_vec[6], eps_vec[0], eps_vec[2], eps_vec[3], eps_vec[13], eps_vec[10], dummy);
}
if (threephase || oilGas) {
// ### krg
initEPSParam(cell, eps_transf[cell].gas, false,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
funcForCell(cell).smax_[gpos],
funcForCell(cell).sogcr_,
oilGas ? -1.0 : funcForCell(cell).smin_[wpos],
funcForCell(cell).krgr_,
funcForCell(cell).krgmax_,
funcForCell(cell).pcgmax_,
satFuncBase.smin_[gpos],
satFuncBase.sgcr_,
satFuncBase.smax_[gpos],
satFuncBase.sogcr_,
oilGas ? -1.0 : satFuncBase.smin_[wpos],
satFuncBase.krgr_,
satFuncBase.krgmax_,
satFuncBase.pcgmax_,
eps_vec[3], eps_vec[5], eps_vec[4], eps_vec[7], eps_vec[0], eps_vec[12], eps_vec[9], eps_vec[16]);
// ### krog
initEPSParam(cell, eps_transf[cell].gasoil, true,
0.0,
funcForCell(cell).sogcr_,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
oilGas ? -1.0 : funcForCell(cell).smin_[wpos],
funcForCell(cell).krorg_,
funcForCell(cell).kromax_,
satFuncBase.sogcr_,
satFuncBase.smin_[gpos],
satFuncBase.sgcr_,
oilGas ? -1.0 : satFuncBase.smin_[wpos],
satFuncBase.krorg_,
satFuncBase.kromax_,
0.0,
eps_vec[3], eps_vec[7], eps_vec[3], eps_vec[5], eps_vec[0], eps_vec[14], eps_vec[10], dummy);
}
@@ -542,9 +567,8 @@ namespace Opm
}
// Initialize saturation scaling parameter
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(Opm::DeckConstPtr deck,
void SaturationPropsFromDeck::initEPSKey(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
@@ -553,9 +577,9 @@ namespace Opm
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];
const bool useAqua = phase_usage_.phase_used[BlackoilPhases::Aqua];
const bool useLiquid = phase_usage_.phase_used[BlackoilPhases::Liquid];
const bool useVapour = phase_usage_.phase_used[BlackoilPhases::Vapour];
bool useKeyword = deck->hasKeyword(keyword);
bool useStateKeyword = eclipseState->hasDoubleGridProperty(keyword);
const std::map<std::string, int> kw2tab = {
@@ -594,49 +618,49 @@ namespace Opm
itab = 1;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwmax_;
scaleparam[i] = satfunc_[cell_to_func_[i]].getSatFuncBase().krwmax_;
}
} else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 1))) {
itab = 2;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgmax_;
scaleparam[i] = satfunc_[cell_to_func_[i]].getSatFuncBase().krgmax_;
}
} else if (keyword == std::string("KRO") || keyword == std::string("IKRO") ) {
if (useLiquid && (useKeyword || columnIsMasked_(deck, "ENKRVD", 2))) {
itab = 3;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).kromax_;
scaleparam[i] = satfunc_[cell_to_func_[i]].getSatFuncBase().kromax_;
}
} else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENKRVD", 3))) {
itab = 4;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwr_;
scaleparam[i] = satfunc_[cell_to_func_[i]].getSatFuncBase().krwr_;
}
} else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 4))) {
itab = 5;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgr_;
scaleparam[i] = satfunc_[cell_to_func_[i]].getSatFuncBase().krgr_;
}
} else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENKRVD", 5))) {
itab = 6;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorw_;
scaleparam[i] = satfunc_[cell_to_func_[i]].getSatFuncBase().krorw_;
}
} else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 6))) {
itab = 7;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorg_;
scaleparam[i] = satfunc_[cell_to_func_[i]].getSatFuncBase().krorg_;
}
} else {
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
@@ -657,11 +681,11 @@ namespace Opm
if (useAqua && (keyword == std::string("PCW") || keyword == std::string("IPCW")) ) {
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).pcwmax_;
scaleparam[i] = satfunc_[cell_to_func_[i]].getSatFuncBase().pcwmax_;
} else if (useVapour && (keyword == std::string("PCG") || keyword == std::string("IPCG")) ) {
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).pcgmax_;
scaleparam[i] = satfunc_[cell_to_func_[i]].getSatFuncBase().pcgmax_;
}
}
@@ -738,8 +762,8 @@ namespace Opm
}
// Saturation scaling
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPSParam(const int cell,
inline
void SaturationPropsFromDeck::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)

94
opm/core/simulator/BlackoilStateToFluidState.hpp Normal file
View File

@@ -0,0 +1,94 @@
/*
Copyright 2015 Andreas Lauser
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOIL_STATE_TO_FLUID_STATE_HEADER_INCLUDED
#define OPM_BLACKOIL_STATE_TO_FLUID_STATE_HEADER_INCLUDED
#include <opm/core/simulator/BlackoilState.hpp>
namespace Opm
{
/*!
* \brief This is an light weight "impedance adaption" class with a well defined API for
* saturation and PVT functions.
*
* It uses a stripped down version of opm-material's FluidState API and takes an
* Opm::BlackoilState plus a cell index.
*
* Note that this class requires that is underlying BlackoilState must valid for at least
* as long as an object of BlackoilStateToFluidState is used.
*/
class BlackoilStateToFluidState
{
public:
typedef double Scalar;
enum { numPhases = 3 };
enum { numComponents = 3 };
/*!
* \brief Create a BlackoilState to Fluid state wrapper object.
*
* Note that this class requires that is underlying BlackoilState must valid for at least
* as long as an object of BlackoilStateToFluidState is used.
*/
BlackoilStateToFluidState(const BlackoilState& blackoilState)
: blackoilState_(blackoilState)
{
if (blackoilState_.numPhases() != numPhases) {
OPM_THROW(std::runtime_error,
"Only " << numPhases << " are supported, but the deck specifies " << blackoilState_.numPhases());
}
}
/*!
* \brief Sets the index of the currently used cell.
*
* After calling this, the values returned by the other methods are specific for this
* cell.
*/
void setCurrentCellIndex(unsigned cellIdx)
{ cellIdx_ = cellIdx; }
/*!
* \brief Returns the saturation of a phase for the current cell index.
*/
Scalar saturation(int phaseIdx) const
{ return blackoilState_.saturation()[numPhases*cellIdx_ + phaseIdx]; }
/*!
* \brief Returns the temperature [K] of a phase for the current cell index.
*/
Scalar temperature(int phaseIdx) const
{ return blackoilState_.temperature()[cellIdx_]; }
// TODO (?) pressure, composition, etc
private:
size_t numCells() const
{ return blackoilState_.pressure().size(); }
const BlackoilState& blackoilState_;
unsigned cellIdx_;
}
} // namespace Opm
#endif // OPM_SIMULATORTIMER_HEADER_INCLUDED

97
opm/core/simulator/ExplicitArraysFluidState.hpp Normal file
View File

@@ -0,0 +1,97 @@
/*
Copyright 2015 Andreas Lauser
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_EXPLICIT_ARRAYS_FLUID_STATE_HEADER_INCLUDED
#define OPM_EXPLICIT_ARRAYS_FLUID_STATE_HEADER_INCLUDED
#include <opm/core/simulator/BlackoilState.hpp>
namespace Opm
{
/*!
* \brief This is a fluid state which translates global arrays and translates them to a
* subset of the fluid state API.
*
* This class is similar to Opm::BlackoilStateToFluidState.
*/
class ExplicitArraysFluidState
{
public:
typedef double Scalar;
enum { numPhases = 3 };
enum { numComponents = 3 };
ExplicitArraysFluidState()
{}
/*!
* \brief Sets the currently used array index.
*
* After calling this, the values returned by the other methods are specific for this
* index.
*/
void setIndex(unsigned arrayIdx)
{ arrayIdx_ = arrayIdx; }
/*!
* \brief Set the array containing the phase saturations.
*
* This array is supposed to be of size numPhases*size and is not allowed to be
* deleted while the ExplicitArraysFluidState object is alive. This class assumes
* that the saturations of all phase saturations for a point are consequtive, i.e.,
* in the array the saturations cycle fastest.
*/
void setSaturationArray(const double* saturations)
{ saturations_ = saturations; }
/*!
* \brief Set the array containing the phase temperatures.
*
* This array is supposed to be of size 'size' and is not allowed to be
* deleted while the ExplicitArraysFluidState object is alive.
*/
void setTemperatureArray(const double* temperature)
{ temperature_ = temperature; }
/*!
* \brief Returns the saturation of a phase for the current cell index.
*/
Scalar saturation(int phaseIdx) const
{ return saturations_[numPhases*arrayIdx_ + phaseIdx]; }
/*!
* \brief Returns the temperature [K] of a phase for the current cell index.
*/
Scalar temperature(int phaseIdx) const
{ return temperature_[arrayIdx_]; }
// TODO (?) pressure, composition, etc
private:
const double* saturations_;
const double* temperature_;
unsigned arrayIdx_;
};
} // namespace Opm
#endif // OPM_SIMULATORTIMER_HEADER_INCLUDED