Merge pull request #834 from andlaus/satfunc_refactoring
use opm-material for the saturation functions
This commit is contained in:
Atgeirr Flø Rasmussen
@@ -98,8 +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/SatFuncSimple.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
|
||||
@@ -363,14 +363,8 @@ list (APPEND PUBLIC_HEADER_FILES
|
||||
opm/core/props/rock/RockBasic.hpp
|
||||
opm/core/props/rock/RockCompressibility.hpp
|
||||
opm/core/props/rock/RockFromDeck.hpp
|
||||
opm/core/props/satfunc/SatFuncBase.hpp
|
||||
opm/core/props/satfunc/SatFuncGwseg.hpp
|
||||
opm/core/props/satfunc/SatFuncMultiplexer.hpp
|
||||
opm/core/props/satfunc/SatFuncSimple.hpp
|
||||
opm/core/props/satfunc/SatFuncStone2.hpp
|
||||
opm/core/props/satfunc/SaturationPropsBasic.hpp
|
||||
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
|
||||
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
|
||||
opm/core/props/satfunc/SaturationPropsInterface.hpp
|
||||
opm/core/simulator/AdaptiveSimulatorTimer.hpp
|
||||
opm/core/simulator/AdaptiveTimeStepping.hpp
|
||||
@@ -379,6 +373,7 @@ list (APPEND PUBLIC_HEADER_FILES
|
||||
opm/core/simulator/BlackoilStateToFluidState.hpp
|
||||
opm/core/simulator/EquilibrationHelpers.hpp
|
||||
opm/core/simulator/ExplicitArraysFluidState.hpp
|
||||
opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp
|
||||
opm/core/simulator/TimeStepControl.hpp
|
||||
opm/core/simulator/SimulatorCompressibleTwophase.hpp
|
||||
opm/core/simulator/SimulatorIncompTwophase.hpp
|
||||
|
||||
4
debian/control
vendored
4
debian/control
vendored
@@ -6,8 +6,8 @@ Build-Depends: build-essential, debhelper (>= 9), libboost-filesystem-dev,
|
||||
libsuperlu3-dev, gfortran, libsuitesparse-dev, pkg-config,
|
||||
libdune-common-dev, libdune-istl-dev, cmake, libtinyxml-dev, bc,
|
||||
libert.ecl-dev, git, zlib1g-dev, libtool, doxygen, libopm-parser-dev,
|
||||
texlive-latex-extra, texlive-latex-recommended, ghostscript,
|
||||
libboost-iostreams-dev
|
||||
libopm-material-dev, texlive-latex-extra, texlive-latex-recommended,
|
||||
ghostscript, libboost-iostreams-dev
|
||||
Standards-Version: 3.9.2
|
||||
Section: libs
|
||||
Homepage: http://opm-project.org
|
||||
|
||||
@@ -4,4 +4,4 @@ Description: Open Porous Media Initiative Core Library
|
||||
Version: 1.1
|
||||
Label: 2013.10
|
||||
Maintainer: atgeirr@sintef.no
|
||||
Depends: dune-common (>= 2.2) dune-istl (>= 2.2) opm-cmake opm-parser
|
||||
Depends: dune-common (>= 2.2) dune-istl (>= 2.2) opm-cmake opm-parser opm-material
|
||||
|
||||
@@ -28,7 +28,20 @@ namespace Opm
|
||||
const UnstructuredGrid& grid,
|
||||
bool init_rock)
|
||||
{
|
||||
init(deck, eclState, grid.number_of_cells, grid.global_cell, grid.cartdims,
|
||||
std::vector<int> compressedToCartesianIdx(grid.number_of_cells);
|
||||
for (unsigned cellIdx = 0; cellIdx < grid.number_of_cells; ++cellIdx) {
|
||||
if (grid.global_cell) {
|
||||
compressedToCartesianIdx[cellIdx] = grid.global_cell[cellIdx];
|
||||
}
|
||||
else {
|
||||
compressedToCartesianIdx[cellIdx] = cellIdx;
|
||||
}
|
||||
}
|
||||
|
||||
auto materialLawManager = std::make_shared<MaterialLawManager>();
|
||||
materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx);
|
||||
|
||||
init(deck, eclState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cartdims,
|
||||
grid.cell_centroids, grid.dimensions, init_rock);
|
||||
}
|
||||
|
||||
@@ -38,7 +51,20 @@ namespace Opm
|
||||
const parameter::ParameterGroup& param,
|
||||
bool init_rock)
|
||||
{
|
||||
init(deck, eclState, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids,
|
||||
std::vector<int> compressedToCartesianIdx(grid.number_of_cells);
|
||||
for (unsigned cellIdx = 0; cellIdx < grid.number_of_cells; ++cellIdx) {
|
||||
if (grid.global_cell) {
|
||||
compressedToCartesianIdx[cellIdx] = grid.global_cell[cellIdx];
|
||||
}
|
||||
else {
|
||||
compressedToCartesianIdx[cellIdx] = cellIdx;
|
||||
}
|
||||
}
|
||||
|
||||
auto materialLawManager = std::make_shared<MaterialLawManager>();
|
||||
materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx);
|
||||
|
||||
init(deck, eclState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids,
|
||||
grid.dimensions, param, init_rock);
|
||||
}
|
||||
|
||||
|
||||
@@ -27,6 +27,8 @@
|
||||
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
|
||||
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
|
||||
#include <memory>
|
||||
@@ -41,6 +43,8 @@ namespace Opm
|
||||
class BlackoilPropertiesFromDeck : public BlackoilPropertiesInterface
|
||||
{
|
||||
public:
|
||||
typedef typename SaturationPropsFromDeck::MaterialLawManager MaterialLawManager;
|
||||
|
||||
/// Initialize from deck and grid.
|
||||
/// \param[in] deck Deck input parser
|
||||
/// \param[in] grid Grid to which property object applies, needed for the
|
||||
@@ -67,7 +71,6 @@ namespace Opm
|
||||
const parameter::ParameterGroup& param,
|
||||
bool init_rock=true);
|
||||
|
||||
|
||||
template<class CentroidIterator>
|
||||
BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclState,
|
||||
@@ -89,6 +92,18 @@ namespace Opm
|
||||
const parameter::ParameterGroup& param,
|
||||
bool init_rock=true);
|
||||
|
||||
template<class CentroidIterator>
|
||||
BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclState,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const int* cart_dims,
|
||||
const CentroidIterator& begin_cell_centroids,
|
||||
int dimension,
|
||||
const parameter::ParameterGroup& param,
|
||||
bool init_rock=true);
|
||||
|
||||
/// Destructor.
|
||||
virtual ~BlackoilPropertiesFromDeck();
|
||||
|
||||
@@ -240,6 +255,7 @@ namespace Opm
|
||||
template<class CentroidIterator>
|
||||
void init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclState,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const int* cart_dims,
|
||||
@@ -249,6 +265,7 @@ namespace Opm
|
||||
template<class CentroidIterator>
|
||||
void init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclState,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const int* cart_dims,
|
||||
@@ -256,9 +273,11 @@ namespace Opm
|
||||
int dimension,
|
||||
const parameter::ParameterGroup& param,
|
||||
bool init_rock);
|
||||
|
||||
RockFromDeck rock_;
|
||||
std::vector<int> cellPvtRegionIdx_;
|
||||
BlackoilPvtProperties pvt_;
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager_;
|
||||
std::shared_ptr<SaturationPropsInterface> satprops_;
|
||||
mutable std::vector<double> B_;
|
||||
mutable std::vector<double> dB_;
|
||||
|
||||
@@ -1,3 +1,5 @@
|
||||
#include <memory>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
template<class CentroidIterator>
|
||||
@@ -10,7 +12,19 @@ namespace Opm
|
||||
int dimension,
|
||||
bool init_rock)
|
||||
{
|
||||
init(deck, eclState, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
|
||||
std::vector<int> compressedToCartesianIdx(number_of_cells);
|
||||
for (unsigned cellIdx = 0; cellIdx < number_of_cells; ++cellIdx) {
|
||||
if (global_cell) {
|
||||
compressedToCartesianIdx[cellIdx] = global_cell[cellIdx];
|
||||
}
|
||||
else {
|
||||
compressedToCartesianIdx[cellIdx] = cellIdx;
|
||||
}
|
||||
}
|
||||
auto materialLawManager = std::make_shared<MaterialLawManager>();
|
||||
materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx);
|
||||
|
||||
init(deck, eclState, materialLawManager, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
|
||||
init_rock);
|
||||
}
|
||||
|
||||
@@ -25,8 +39,45 @@ namespace Opm
|
||||
const parameter::ParameterGroup& param,
|
||||
bool init_rock)
|
||||
{
|
||||
std::vector<int> compressedToCartesianIdx(number_of_cells);
|
||||
for (unsigned cellIdx = 0; cellIdx < number_of_cells; ++cellIdx) {
|
||||
if (global_cell) {
|
||||
compressedToCartesianIdx[cellIdx] = global_cell[cellIdx];
|
||||
}
|
||||
else {
|
||||
compressedToCartesianIdx[cellIdx] = cellIdx;
|
||||
}
|
||||
}
|
||||
auto materialLawManager = std::make_shared<MaterialLawManager>();
|
||||
materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx);
|
||||
|
||||
init(deck,
|
||||
eclState,
|
||||
materialLawManager,
|
||||
number_of_cells,
|
||||
global_cell,
|
||||
cart_dims,
|
||||
begin_cell_centroids,
|
||||
dimension,
|
||||
param,
|
||||
init_rock);
|
||||
}
|
||||
|
||||
template<class CentroidIterator>
|
||||
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclState,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const int* cart_dims,
|
||||
const CentroidIterator& begin_cell_centroids,
|
||||
int dimension,
|
||||
const parameter::ParameterGroup& param,
|
||||
bool init_rock)
|
||||
{
|
||||
init(deck,
|
||||
eclState,
|
||||
materialLawManager,
|
||||
number_of_cells,
|
||||
global_cell,
|
||||
cart_dims,
|
||||
@@ -39,6 +90,7 @@ namespace Opm
|
||||
template<class CentroidIterator>
|
||||
inline void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclState,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const int* cart_dims,
|
||||
@@ -56,8 +108,8 @@ namespace Opm
|
||||
pvt_.init(deck, eclState, /*numSamples=*/0);
|
||||
SaturationPropsFromDeck* ptr
|
||||
= new SaturationPropsFromDeck();
|
||||
ptr->init(phaseUsageFromDeck(deck), materialLawManager);
|
||||
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 ("
|
||||
@@ -68,6 +120,7 @@ namespace Opm
|
||||
template<class CentroidIterator>
|
||||
inline void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclState,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const int* cart_dims,
|
||||
@@ -95,8 +148,8 @@ namespace Opm
|
||||
|
||||
SaturationPropsFromDeck* ptr
|
||||
= new SaturationPropsFromDeck();
|
||||
ptr->init(phaseUsageFromDeck(deck), materialLawManager);
|
||||
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 ("
|
||||
|
||||
@@ -32,7 +32,20 @@ namespace Opm
|
||||
{
|
||||
rock_.init(eclState, grid.number_of_cells, grid.global_cell, grid.cartdims);
|
||||
pvt_.init(deck);
|
||||
satprops_.init(deck, eclState, grid);
|
||||
auto materialLawManager = std::make_shared<typename SaturationPropsFromDeck::MaterialLawManager>();
|
||||
|
||||
std::vector<int> compressedToCartesianIdx(grid.number_of_cells);
|
||||
for (unsigned cellIdx = 0; cellIdx < grid.number_of_cells; ++cellIdx) {
|
||||
if (grid.global_cell) {
|
||||
compressedToCartesianIdx[cellIdx] = grid.global_cell[cellIdx];
|
||||
}
|
||||
else {
|
||||
compressedToCartesianIdx[cellIdx] = cellIdx;
|
||||
}
|
||||
}
|
||||
materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx);
|
||||
|
||||
satprops_.init(deck, eclState, materialLawManager, 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() << ").");
|
||||
|
||||
@@ -1,329 +0,0 @@
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#ifndef SATFUNCBASE_HPP
|
||||
#define SATFUNCBASE_HPP
|
||||
|
||||
#include <opm/core/utility/UniformTableLinear.hpp>
|
||||
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
|
||||
#include <opm/core/utility/NonuniformTableLinear.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
// Transforms for saturation table scaling
|
||||
struct EPSTransforms {
|
||||
struct Transform {
|
||||
bool doNotScale;
|
||||
bool do_3pt;
|
||||
double smin;
|
||||
double scr;
|
||||
double sr;
|
||||
double smax;
|
||||
double slope1;
|
||||
double slope2;
|
||||
double scaleSat(double ss, double s_r, double s_cr, double s_max) const;
|
||||
double scaleSatInv(double s, double s_r, double s_cr, double s_max) const;
|
||||
double scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const; // Returns scaleSat'(s)
|
||||
double scaleSatPc(double s, double s_min, double s_max) const;
|
||||
double scaleSatDerivPc(double s, double s_min, double s_max) const; // Returns scaleSatPc'(s)
|
||||
bool doKrMax;
|
||||
bool doKrCrit;
|
||||
bool doSatInterp;
|
||||
double krsr;
|
||||
double krmax;
|
||||
double krSlopeMax;
|
||||
double krSlopeCrit;
|
||||
double scaleKr(double s, double kr, double krsr_tab) const;
|
||||
double scaleKrDeriv(double s, double krDeriv) const; // Returns scaleKr'(kr(scaleSat(s)))*kr'((scaleSat(s))
|
||||
double pcFactor; // Scaling factor for capillary pressure.
|
||||
void printMe(std::ostream & out);
|
||||
};
|
||||
|
||||
Transform wat;
|
||||
Transform watoil;
|
||||
Transform gas;
|
||||
Transform gasoil;
|
||||
};
|
||||
|
||||
// Hysteresis
|
||||
struct SatHyst {
|
||||
SatHyst() {
|
||||
sg_hyst = -1.0;
|
||||
sow_hyst = -1.0;
|
||||
|
||||
sg_shift = 0.0;
|
||||
sow_shift = 0.0;
|
||||
}
|
||||
|
||||
double sg_hyst;
|
||||
double sg_shift;
|
||||
double sow_hyst;
|
||||
double sow_shift;
|
||||
void printMe(std::ostream & out);
|
||||
};
|
||||
|
||||
|
||||
template <class TableType>
|
||||
class SatFuncBase : public BlackoilPhases
|
||||
{
|
||||
public:
|
||||
virtual ~SatFuncBase() {};
|
||||
|
||||
void init(Opm::EclipseStateConstPtr eclipseState,
|
||||
const int table_num,
|
||||
const PhaseUsage phase_usg,
|
||||
const int samples);
|
||||
void updateSatHyst(const double* s,
|
||||
const EPSTransforms* epst,
|
||||
const EPSTransforms* epst_hyst,
|
||||
SatHyst* sat_hyst) const;
|
||||
double smin_[PhaseUsage::MaxNumPhases];
|
||||
double smax_[PhaseUsage::MaxNumPhases];
|
||||
double krwmax_; // Max water relperm
|
||||
double krgmax_; // Max gas relperm
|
||||
double kromax_; // Max oil relperm
|
||||
double swcr_; // Critical water saturation.
|
||||
double sgcr_; // Critical gas saturation.
|
||||
double krwr_; // Water relperm at critical oil-in-water saturation.
|
||||
double krgr_; // Gas relperm at critical oil-in-gas saturation.
|
||||
double sowcr_; // Critical oil-in-water saturation.
|
||||
double sogcr_; // Critical oil-in-gas-and-connate-water saturation.
|
||||
double krorw_; // Oil relperm at critical water saturation.
|
||||
double krorg_; // Oil relperm at critical gas saturation.
|
||||
double pcwmax_; // Max oil-water capillary pressure.
|
||||
double pcgmax_; // Max gas-oil capillary pressure.
|
||||
|
||||
protected:
|
||||
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
|
||||
TableType krw_;
|
||||
TableType krow_;
|
||||
TableType pcow_;
|
||||
TableType krg_;
|
||||
TableType krog_; TableType pcog_;
|
||||
double krocw_; // = krow_(s_wc)
|
||||
|
||||
private:
|
||||
void extendTable(const std::vector<double>& xv,
|
||||
std::vector<double>& xv_ex,
|
||||
double pm) const;
|
||||
void initializeTableType(TableType& table,
|
||||
const std::vector<double>& arg,
|
||||
const std::vector<double>& value,
|
||||
const int samples);
|
||||
};
|
||||
|
||||
template <class TableType>
|
||||
void SatFuncBase<TableType>::init(Opm::EclipseStateConstPtr eclipseState,
|
||||
const int table_num,
|
||||
const PhaseUsage phase_usg,
|
||||
const int samples)
|
||||
{
|
||||
std::shared_ptr<const TableManager> tables = eclipseState->getTableManager();
|
||||
phase_usage = phase_usg;
|
||||
double swco = 0.0;
|
||||
double swmax = 1.0;
|
||||
if (phase_usage.phase_used[Aqua]) {
|
||||
const Opm::SwofTable& swof(tables->getSwofTables()[table_num]);
|
||||
const std::vector<double>& sw = swof.getSwColumn();
|
||||
const std::vector<double>& krw = swof.getKrwColumn();
|
||||
const std::vector<double>& krow = swof.getKrowColumn();
|
||||
const std::vector<double>& pcow = swof.getPcowColumn();
|
||||
if (krw.front() != 0.0 || krow.back() != 0.0) {
|
||||
OPM_THROW(std::runtime_error, "Error SWOF data - non-zero krw(swco) and/or krow(1-sor)");
|
||||
}
|
||||
|
||||
// Extend the tables with constant values such that the
|
||||
// derivatives at the endpoints are zero
|
||||
int n = sw.size();
|
||||
std::vector<double> sw_ex(n+2);
|
||||
std::vector<double> krw_ex(n+2);
|
||||
std::vector<double> krow_ex(n+2);
|
||||
std::vector<double> pcow_ex(n+2);
|
||||
|
||||
extendTable(sw,sw_ex,1);
|
||||
extendTable(krw,krw_ex,0);
|
||||
extendTable(krow,krow_ex,0);
|
||||
extendTable(pcow,pcow_ex,0);
|
||||
|
||||
initializeTableType(krw_,sw_ex, krw_ex, samples);
|
||||
initializeTableType(krow_,sw_ex, krow_ex, samples);
|
||||
initializeTableType(pcow_,sw_ex, pcow_ex, samples);
|
||||
|
||||
krocw_ = krow[0]; // At connate water -> ecl. SWOF
|
||||
swco = sw[0];
|
||||
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
|
||||
swmax = sw.back();
|
||||
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
|
||||
|
||||
krwmax_ = krw.back();
|
||||
kromax_ = krow.front();
|
||||
swcr_ = swmax;
|
||||
sowcr_ = 1.0 - swco;
|
||||
krwr_ = krw.back();
|
||||
krorw_ = krow.front();
|
||||
for (std::vector<double>::size_type i=1; i<sw.size(); ++i) {
|
||||
if (krw[i]> 0.0) {
|
||||
swcr_ = sw[i-1];
|
||||
krorw_ = krow[i-1];
|
||||
break;
|
||||
}
|
||||
}
|
||||
for (std::vector<double>::size_type i=sw.size()-1; i>=1; --i) {
|
||||
if (krow[i-1]> 0.0) {
|
||||
sowcr_ = 1.0 - sw[i];
|
||||
krwr_ = krw[i];
|
||||
break;
|
||||
}
|
||||
}
|
||||
pcwmax_ = pcow.front();
|
||||
}
|
||||
if (phase_usage.phase_used[Vapour]) {
|
||||
const Opm::SgofTable& sgof = tables->getSgofTables()[table_num];
|
||||
const std::vector<double>& sg = sgof.getSgColumn();
|
||||
const std::vector<double>& krg = sgof.getKrgColumn();
|
||||
const std::vector<double>& krog = sgof.getKrogColumn();
|
||||
const std::vector<double>& pcog = sgof.getPcogColumn();
|
||||
|
||||
// Extend the tables with constant values such that the
|
||||
// derivatives at the endpoints are zero
|
||||
int n = sg.size();
|
||||
std::vector<double> sg_ex(n+2);
|
||||
std::vector<double> krg_ex(n+2);
|
||||
std::vector<double> krog_ex(n+2);
|
||||
std::vector<double> pcog_ex(n+2);
|
||||
|
||||
extendTable(sg,sg_ex,1);
|
||||
extendTable(krg,krg_ex,0);
|
||||
extendTable(krog,krog_ex,0);
|
||||
extendTable(pcog,pcog_ex,0);
|
||||
|
||||
initializeTableType(krg_,sg_ex, krg_ex, samples);
|
||||
initializeTableType(krog_,sg_ex, krog_ex, samples);
|
||||
initializeTableType(pcog_,sg_ex, pcog_ex, samples);
|
||||
|
||||
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
|
||||
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
|
||||
OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() <<
|
||||
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
|
||||
}
|
||||
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
|
||||
smin_[phase_usage.phase_pos[Vapour]] = sg.front();
|
||||
krgmax_ = krg.back();
|
||||
|
||||
sgcr_ = sg.front();
|
||||
sogcr_ = 1.0 - sg.back();
|
||||
krgr_ = krg.back();
|
||||
krorg_ = krg.front();
|
||||
for (std::vector<double>::size_type i=1; i<sg.size(); ++i) {
|
||||
if (krg[i]> 0.0) {
|
||||
sgcr_ = sg[i-1];
|
||||
krorg_ = krog[i-1];
|
||||
break;
|
||||
}
|
||||
}
|
||||
for (std::vector<double>::size_type i=sg.size()-1; i>=1; --i) {
|
||||
if (krog[i-1]> 0.0) {
|
||||
sogcr_ = 1.0 - sg[i];
|
||||
krgr_ = krg[i];
|
||||
break;
|
||||
}
|
||||
}
|
||||
pcgmax_ = pcog.back();
|
||||
|
||||
}
|
||||
|
||||
if (phase_usage.phase_used[Vapour] && phase_usage.phase_used[Aqua]) {
|
||||
sowcr_ -= smin_[phase_usage.phase_pos[Vapour]];
|
||||
sogcr_ -= smin_[phase_usage.phase_pos[Aqua]];
|
||||
smin_[phase_usage.phase_pos[Liquid]] = 0.0;
|
||||
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]]
|
||||
- smin_[phase_usage.phase_pos[Vapour]]; // First entry in SGOF-table supposed to be zero anyway ...
|
||||
} else if (phase_usage.phase_used[Aqua]) {
|
||||
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Aqua]];
|
||||
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]];
|
||||
} else if (phase_usage.phase_used[Vapour]) {
|
||||
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Vapour]];
|
||||
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Vapour]];
|
||||
}
|
||||
}
|
||||
|
||||
template <class TableType>
|
||||
void SatFuncBase<TableType>::updateSatHyst(const double* s,
|
||||
const EPSTransforms* epst,
|
||||
const EPSTransforms* epst_hyst,
|
||||
SatHyst* sat_hyst) const
|
||||
{
|
||||
if (phase_usage.phase_used[Aqua] && phase_usage.phase_used[Vapour]) { //Water/Oil/Gas
|
||||
int opos = phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
if (s[opos] > sat_hyst->sow_hyst)
|
||||
{
|
||||
sat_hyst->sow_hyst = s[opos];
|
||||
double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]);
|
||||
double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]);
|
||||
sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst;
|
||||
}
|
||||
if (s[gpos] > sat_hyst->sg_hyst)
|
||||
{
|
||||
sat_hyst->sg_hyst = s[gpos];
|
||||
double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]);
|
||||
double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]);
|
||||
sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst;
|
||||
}
|
||||
} else if (phase_usage.phase_used[Aqua]) { //Water/oil
|
||||
int opos = phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
if (s[opos] > sat_hyst->sow_hyst)
|
||||
{
|
||||
sat_hyst->sow_hyst = s[opos];
|
||||
double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]);
|
||||
double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]);
|
||||
sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst;
|
||||
}
|
||||
} else if (phase_usage.phase_used[Vapour]) {//Gas/Oil
|
||||
int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
if (s[gpos] > sat_hyst->sg_hyst)
|
||||
{
|
||||
sat_hyst->sg_hyst = s[gpos];
|
||||
double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]);
|
||||
double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]);
|
||||
sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class TableType>
|
||||
void SatFuncBase<TableType>::extendTable(const std::vector<double>& xv,
|
||||
std::vector<double>& xv_ex,
|
||||
double pm) const
|
||||
{
|
||||
int n = xv.size();
|
||||
xv_ex[0] = xv[0]-pm;
|
||||
xv_ex[n+1] = xv[n-1]+pm;
|
||||
for (int i=0; i<n; i++)
|
||||
{
|
||||
xv_ex[i+1] = xv[i];
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
#endif // SATFUNCBASE_HPP
|
||||
@@ -1,637 +0,0 @@
|
||||
/*
|
||||
Copyright 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#ifndef SATFUNCGWSEG_HPP
|
||||
#define SATFUNCGWSEG_HPP
|
||||
|
||||
#include <opm/core/props/satfunc/SatFuncBase.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
template<class TableType>
|
||||
class SatFuncGwseg : public SatFuncBase<TableType>
|
||||
{
|
||||
public:
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& fluidState, double* kr) const;
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const;
|
||||
template <class FluidState>
|
||||
void evalPc(const FluidState& fluidState, double* pc) const;
|
||||
template <class FluidState>
|
||||
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const;
|
||||
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst) const;
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const;
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const;
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const;
|
||||
template <class FluidState>
|
||||
void evalPc(const FluidState& fluidState, double* pc, const EPSTransforms* epst) const;
|
||||
template <class FluidState>
|
||||
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds, const EPSTransforms* epst) const;
|
||||
};
|
||||
|
||||
typedef SatFuncGwseg<UniformTableLinear<double> > SatFuncGwsegUniform;
|
||||
typedef SatFuncGwseg<NonuniformTableLinear<double> > SatFuncGwsegNonuniform;
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr) const
|
||||
{
|
||||
if (this->phase_usage.num_phases == 3) {
|
||||
// Relative permeability model based on segregation of water
|
||||
// and gas, with oil present in both water and gas zones.
|
||||
double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
|
||||
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
|
||||
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
|
||||
const double eps = 1e-5;
|
||||
swco = std::min(swco,sw-eps);
|
||||
|
||||
// xw and xg are the fractions occupied by water and gas zones.
|
||||
const double ssg = sw - swco + sg;
|
||||
const double xw = (sw - swco) / ssg;
|
||||
const double xg = 1 - xw;
|
||||
const double ssw = sg + sw;
|
||||
|
||||
const double krw = this->krw_(sw);
|
||||
const double krg = this->krg_(sg);
|
||||
const double krow = this->krow_(ssw);
|
||||
const double krog = this->krog_(ssg);
|
||||
kr[BlackoilPhases::Aqua] = krw;
|
||||
kr[BlackoilPhases::Vapour] = krg;
|
||||
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
|
||||
return;
|
||||
}
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double krow = this->krow_(sw);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double krog = this->krog_(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst) const
|
||||
{
|
||||
if (this->phase_usage.num_phases == 3) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
|
||||
// Relative permeability model based on segregation of water
|
||||
// and gas, with oil present in both water and gas zones.
|
||||
|
||||
// TODO Also consider connate gas ...
|
||||
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
|
||||
double swco = epst->wat.smin;
|
||||
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
|
||||
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
|
||||
const double eps = 1e-6;
|
||||
swco = std::min(swco,sw-eps);
|
||||
const double ssw = sg + sw;
|
||||
const double ssg = std::max(sg + sw - swco, eps);
|
||||
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
|
||||
|
||||
double ssow = 1.0-ssw;
|
||||
double ssog = 1.0-ssg-swco;
|
||||
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
|
||||
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
|
||||
const double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
|
||||
const double krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
|
||||
const double krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
|
||||
const double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
|
||||
|
||||
// xw and xg are the fractions occupied by water and gas zones.
|
||||
const double xw = (sw - swco) / d;
|
||||
const double xg = 1 - xw;
|
||||
kr[BlackoilPhases::Aqua] = krw;
|
||||
kr[BlackoilPhases::Vapour] = krg;
|
||||
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
|
||||
return;
|
||||
}
|
||||
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double krow = this->krow_(sw);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double krog = this->krog_(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
}
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncGwseg<TableType>::evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
|
||||
{
|
||||
if (this->phase_usage.num_phases == 3) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
|
||||
// Relative permeability model based on segregation of water
|
||||
// and gas, with oil present in both water and gas zones.
|
||||
|
||||
// TODO Consider connate gas ...
|
||||
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
|
||||
double swco = epst->wat.smin;
|
||||
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
|
||||
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
|
||||
const double eps = 1e-6;
|
||||
swco = std::min(swco,sw-eps);
|
||||
const double ssw = sg + sw;
|
||||
const double ssg = std::max(sg + sw - swco, eps);
|
||||
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
|
||||
|
||||
double ssow = 1.0-ssw;
|
||||
double ssog = 1.0-ssg-swco;
|
||||
|
||||
// The code below corresponds to EHYSTR * 0 * * KR/
|
||||
// - wettability properties water>oil>gas.
|
||||
// - Carlsen hysteresis model for non-wetting (scanning=shifted_imb). No hysteresis for wetting phase.
|
||||
// The imb-curve currently only differs from drainage curves via endpoint scaling ...
|
||||
|
||||
// Water - use drainage curve only
|
||||
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
|
||||
double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
|
||||
|
||||
// Gas
|
||||
double krg;
|
||||
if (sg >= sat_hyst->sg_hyst) { // Drainage
|
||||
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
|
||||
} else { // Imbibition
|
||||
double sg_shifted = sg + sat_hyst->sg_shift;
|
||||
double _sg = epst_hyst->gas.scaleSat(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
krg = epst_hyst->gas.scaleKr(sg_shifted, this->krg_(_sg), this->krgr_);
|
||||
}
|
||||
|
||||
// Oil in water
|
||||
double krow;
|
||||
if (ssow >= sat_hyst->sow_hyst) { // Drainage
|
||||
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
|
||||
} else { // Imbibition
|
||||
double ssow_shifted = ssow + sat_hyst->sow_shift;
|
||||
double _ssow = epst_hyst->watoil.scaleSat(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(1.0-_ssow), this->krorw_);
|
||||
}
|
||||
|
||||
// Oil in gas and connate water - use drainage curve only
|
||||
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
|
||||
|
||||
// xw and xg are the fractions occupied by water and gas zones.
|
||||
const double xw = (sw - swco) / d;
|
||||
const double xg = 1 - xw;
|
||||
|
||||
// relperms
|
||||
kr[BlackoilPhases::Aqua] = krw;
|
||||
kr[BlackoilPhases::Vapour] = krg;
|
||||
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
|
||||
|
||||
return;
|
||||
}
|
||||
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double krow = this->krow_(sw);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double krog = this->krog_(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
}
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template<class FluidState>
|
||||
void SatFuncGwseg<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
|
||||
{
|
||||
const int np = this->phase_usage.num_phases;
|
||||
std::fill(dkrds, dkrds + np*np, 0.0);
|
||||
|
||||
if (np == 3) {
|
||||
// Relative permeability model based on segregation of water
|
||||
// and gas, with oil present in both water and gas zones.
|
||||
|
||||
double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
|
||||
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
|
||||
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
|
||||
const double eps = 1e-5;
|
||||
swco = std::min(swco,sw-eps);
|
||||
|
||||
// xw and xg are the fractions occupied by water and gas zones.
|
||||
const double ssw = sg + sw;
|
||||
// d = ssg = sw - swco + sg (using 'd' for consistency with mrst docs).
|
||||
const double d = sg + sw - swco;
|
||||
const double xw = (sw - swco) / d;
|
||||
|
||||
const double krw = this->krw_(sw);
|
||||
const double krg = this->krg_(sg);
|
||||
const double krow = this->krow_(ssw);
|
||||
const double krog = this->krog_(d);
|
||||
|
||||
const double xg = 1 - xw;
|
||||
kr[BlackoilPhases::Aqua] = krw;
|
||||
kr[BlackoilPhases::Vapour] = krg;
|
||||
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
|
||||
|
||||
// Derivatives.
|
||||
const double dkrww = this->krw_.derivative(sw);
|
||||
const double dkrgg = this->krg_.derivative(sg);
|
||||
const double dkrow = this->krow_.derivative(ssw);
|
||||
const double dkrog = this->krog_.derivative(d);
|
||||
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
|
||||
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
|
||||
return;
|
||||
}
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double dkrww = this->krw_.derivative(sw);
|
||||
double krow = this->krow_(sw);
|
||||
double dkrow = this->krow_.derivative(sw);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
dkrds[wpos + wpos*np] = dkrww;
|
||||
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double dkrgg = this->krg_.derivative(sg);
|
||||
double krog = this->krog_(sg);
|
||||
double dkrog = this->krog_.derivative(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
dkrds[gpos + gpos*np] = dkrgg;
|
||||
dkrds[opos + gpos*np] = dkrog;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template<class FluidState>
|
||||
void SatFuncGwseg<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const
|
||||
{
|
||||
const int np = this->phase_usage.num_phases;
|
||||
std::fill(dkrds, dkrds + np*np, 0.0);
|
||||
|
||||
if (np == 3) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
|
||||
// Relative permeability model based on segregation of water
|
||||
// and gas, with oil present in both water and gas zones.
|
||||
|
||||
// TODO Also consider connate gas ...
|
||||
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
|
||||
double swco = epst->wat.smin;
|
||||
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
|
||||
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
|
||||
const double eps = 1e-6;
|
||||
swco = std::min(swco,sw-eps);
|
||||
const double ssw = sg + sw;
|
||||
const double ssg = std::max(sg + sw - swco, eps);
|
||||
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
|
||||
|
||||
double ssow = 1.0-ssw;
|
||||
double ssog = 1.0-ssg-swco;
|
||||
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
|
||||
double _dsdsw = epst->wat.scaleSatDeriv(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
|
||||
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
double _dsdsg = epst->gas.scaleSatDeriv(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
double _dsdssow = epst->watoil.scaleSatDeriv(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
double _dsdssog = epst->gasoil.scaleSatDeriv(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
|
||||
const double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
|
||||
const double krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
|
||||
const double krow = epst->watoil.scaleKr(ssow, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
|
||||
const double krog = epst->gasoil.scaleKr(ssog, this->krog_(std::max(1.0-_ssog-_swco,eps)), this->krorg_);
|
||||
|
||||
// xw and xg are the fractions occupied by water and gas zones.
|
||||
const double xw = (sw - swco) / d;
|
||||
const double xg = 1 - xw;
|
||||
kr[BlackoilPhases::Aqua] = krw;
|
||||
kr[BlackoilPhases::Vapour] = krg;
|
||||
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
|
||||
|
||||
// Derivatives.
|
||||
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(sw, this->krw_.derivative(_sw));
|
||||
double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(sg, this->krg_.derivative(_sg));
|
||||
double dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
|
||||
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(std::max(1.0-_ssog-_swco,eps)));
|
||||
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
|
||||
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
|
||||
return;
|
||||
}
|
||||
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double dkrww = this->krw_.derivative(sw);
|
||||
double krow = this->krow_(sw);
|
||||
double dkrow = this->krow_.derivative(sw);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
dkrds[wpos + wpos*np] = dkrww;
|
||||
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double dkrgg = this->krg_.derivative(sg);
|
||||
double krog = this->krog_(sg);
|
||||
double dkrog = this->krog_.derivative(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
dkrds[gpos + gpos*np] = dkrgg;
|
||||
dkrds[opos + gpos*np] = dkrog;
|
||||
}
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template<class FluidState>
|
||||
void SatFuncGwseg<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
|
||||
{
|
||||
const int np = this->phase_usage.num_phases;
|
||||
std::fill(dkrds, dkrds + np*np, 0.0);
|
||||
|
||||
if (np == 3) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
|
||||
// Relative permeability model based on segregation of water
|
||||
// and gas, with oil present in both water and gas zones.
|
||||
|
||||
// TODO Also consider connate gas ...
|
||||
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
|
||||
double swco = epst->wat.smin;
|
||||
const double sw = std::max(fluidState.saturation(BlackoilPhases::Aqua), swco);
|
||||
const double sg = fluidState.saturation(BlackoilPhases::Vapour);
|
||||
const double eps = 1e-6;
|
||||
swco = std::min(swco,sw-eps);
|
||||
const double ssw = sg + sw;
|
||||
const double ssg = std::max(sg + sw - swco, eps);
|
||||
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
|
||||
|
||||
double ssow = 1.0-ssw;
|
||||
double ssog = 1.0-ssg-swco;
|
||||
|
||||
// The code below corresponds to EHYSTR * 0 * * KR/
|
||||
// - wettability properties water>oil>gas.
|
||||
// - Carlsen hysteresis model for non-wetting (scanning=shifted_imb). No hysteresis for wetting phase.
|
||||
// The imb-curve currently only differs from drainage curves via endpoint scaling ...
|
||||
|
||||
// Water - use drainage curve only
|
||||
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
|
||||
double _dsdsw = epst->wat.scaleSatDeriv(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
|
||||
double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
|
||||
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(sw, this->krw_.derivative(_sw));
|
||||
|
||||
// Gas
|
||||
double krg, dkrgg;
|
||||
if (sg >= sat_hyst->sg_hyst) { // Drainage
|
||||
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
double _dsdsg = epst->gas.scaleSatDeriv(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
|
||||
dkrgg = _dsdsg*epst->gas.scaleKrDeriv(sg, this->krg_.derivative(_sg));
|
||||
} else { // Imbibition
|
||||
double sg_shifted = sg + sat_hyst->sg_shift;
|
||||
double _sg = epst_hyst->gas.scaleSat(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
double _dsdsg = epst_hyst->gas.scaleSatDeriv(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
krg = epst_hyst->gas.scaleKr(sg_shifted, this->krg_(_sg), this->krgr_);
|
||||
dkrgg = _dsdsg*epst_hyst->gas.scaleKrDeriv(sg_shifted, this->krg_.derivative(_sg));
|
||||
}
|
||||
|
||||
// Oil in water
|
||||
double krow, dkrow;
|
||||
if (ssow >= sat_hyst->sow_hyst) { // Drainage
|
||||
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
double _dsdssow = epst->watoil.scaleSatDeriv(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
krow = epst->watoil.scaleKr(ssow, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
|
||||
dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
|
||||
} else { // Imbibition
|
||||
double ssow_shifted = ssow + sat_hyst->sow_shift;
|
||||
double _ssow = epst_hyst->watoil.scaleSat(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
double _dsdssow = epst_hyst->watoil.scaleSatDeriv(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
|
||||
dkrow = _dsdssow*epst_hyst->watoil.scaleKrDeriv(ssow_shifted, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
|
||||
}
|
||||
|
||||
// Oil in gas and connate water - use drainage curve only
|
||||
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
double _dsdssog = epst->gasoil.scaleSatDeriv(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
double krog = epst->gasoil.scaleKr(ssog, this->krog_(std::max(1.0-_ssog-_swco,eps)), this->krorg_);
|
||||
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(std::max(1.0-_ssog-_swco,eps)));
|
||||
|
||||
// xw and xg are the fractions occupied by water and gas zones.
|
||||
const double xw = (sw - swco) / d;
|
||||
const double xg = 1 - xw;
|
||||
|
||||
// relperms
|
||||
kr[BlackoilPhases::Aqua] = krw;
|
||||
kr[BlackoilPhases::Vapour] = krg;
|
||||
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
|
||||
|
||||
// Derivatives.
|
||||
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
|
||||
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
|
||||
return;
|
||||
}
|
||||
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double dkrww = this->krw_.derivative(sw);
|
||||
double krow = this->krow_(sw);
|
||||
double dkrow = this->krow_.derivative(sw);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
dkrds[wpos + wpos*np] = dkrww;
|
||||
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double dkrgg = this->krg_.derivative(sg);
|
||||
double krog = this->krog_(sg);
|
||||
double dkrog = this->krog_.derivative(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
dkrds[gpos + gpos*np] = dkrgg;
|
||||
dkrds[opos + gpos*np] = dkrog;
|
||||
}
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncGwseg<TableType>::evalPc(const FluidState& fluidState, double* pc) const
|
||||
{
|
||||
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
pc[pos] = this->pcow_(fluidState.saturation(pos));
|
||||
}
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
pc[pos] = this->pcog_(fluidState.saturation(pos));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncGwseg<TableType>::evalPc(const FluidState& fluidState, double* pc, const EPSTransforms* epst) const
|
||||
{
|
||||
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
double _sw = epst->wat.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
|
||||
pc[pos] = epst->wat.pcFactor*this->pcow_(_sw);
|
||||
}
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
double _sg = epst->gas.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
|
||||
pc[pos] = epst->gas.pcFactor*this->pcog_(_sg);
|
||||
}
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncGwseg<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
|
||||
{
|
||||
// The problem of determining three-phase capillary pressures
|
||||
// is very hard experimentally, usually one extends two-phase
|
||||
// data (as for relative permeability).
|
||||
// In our approach the derivative matrix is quite sparse, only
|
||||
// the diagonal elements corresponding to non-oil phases are
|
||||
// (potentially) nonzero.
|
||||
const int np = this->phase_usage.num_phases;
|
||||
std::fill(dpcds, dpcds + np*np, 0.0);
|
||||
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
pc[pos] = this->pcow_(fluidState.saturation(pos));
|
||||
dpcds[np*pos + pos] = this->pcow_.derivative(fluidState.saturation(pos));
|
||||
}
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
pc[pos] = this->pcog_(fluidState.saturation(pos));
|
||||
dpcds[np*pos + pos] = this->pcog_.derivative(fluidState.saturation(pos));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncGwseg<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds, const EPSTransforms* epst) const
|
||||
{
|
||||
// The problem of determining three-phase capillary pressures
|
||||
// is very hard experimentally, usually one extends two-phase
|
||||
// data (as for relative permeability).
|
||||
// In our approach the derivative matrix is quite sparse, only
|
||||
// the diagonal elements corresponding to non-oil phases are
|
||||
// (potentially) nonzero.
|
||||
const int np = this->phase_usage.num_phases;
|
||||
std::fill(dpcds, dpcds + np*np, 0.0);
|
||||
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
double _sw = epst->wat.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
|
||||
pc[pos] = epst->wat.pcFactor*this->pcow_(_sw);
|
||||
double _dsdsw = epst->wat.scaleSatDerivPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
|
||||
dpcds[np*pos + pos] = epst->wat.pcFactor*_dsdsw*this->pcow_.derivative(_sw);
|
||||
}
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
double _sg = epst->gas.scaleSatPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
|
||||
pc[pos] = epst->gas.pcFactor*this->pcog_(_sg);
|
||||
double _dsdsg = epst->gas.scaleSatDerivPc(fluidState.saturation(pos), this->smin_[pos], this->smax_[pos]);
|
||||
dpcds[np*pos + pos] = epst->gas.pcFactor*_dsdsg*this->pcog_.derivative(_sg);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
#endif // SATFUNCGWSEG_HPP
|
||||
@@ -1,268 +0,0 @@
|
||||
/*
|
||||
Copyright 2015 Andreas Lauser
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
* \copydoc Opm::SatFuncMultiplexer
|
||||
*/
|
||||
#ifndef OPM_SAT_FUNC_MULTIPLEXER_HPP
|
||||
#define OPM_SAT_FUNC_MULTIPLEXER_HPP
|
||||
|
||||
#include "SatFuncGwseg.hpp"
|
||||
#include "SatFuncSimple.hpp"
|
||||
#include "SatFuncStone2.hpp"
|
||||
|
||||
#include <opm/core/props/phaseUsageFromDeck.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
||||
|
||||
#include <memory>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
/*!
|
||||
* \brief This is a multiplexer class which allows to saturation function to be used at
|
||||
* runtime.
|
||||
*/
|
||||
class SatFuncMultiplexer
|
||||
{
|
||||
public:
|
||||
typedef NonuniformTableLinear<double> TableType;
|
||||
|
||||
enum SatFuncType {
|
||||
Gwseg, // <- supposed to be the default
|
||||
Stone2,
|
||||
Simple
|
||||
};
|
||||
|
||||
// this is a helper macro which helps to save us from RSI in the following
|
||||
#define OPM_MULTIPLEXER_CONST const
|
||||
#define OPM_MULTIPLEXER_NON_CONST
|
||||
#define OPM_MULTIPLEXER_SATFUNC_CALL(codeToCall, CONST) \
|
||||
switch (satFuncType_) { \
|
||||
case Gwseg: { \
|
||||
typedef SatFuncGwseg<TableType> SatFunc; \
|
||||
__attribute__((unused)) CONST SatFunc& satFunc = \
|
||||
*static_cast<SatFunc*>(satFunc_.get()); \
|
||||
codeToCall; \
|
||||
break; \
|
||||
} \
|
||||
case Stone2: { \
|
||||
typedef SatFuncStone2<TableType> SatFunc; \
|
||||
__attribute__((unused)) CONST SatFunc& satFunc = \
|
||||
*static_cast<SatFunc*>(satFunc_.get()); \
|
||||
codeToCall; \
|
||||
break; \
|
||||
} \
|
||||
case Simple: { \
|
||||
typedef SatFuncSimple<TableType> SatFunc; \
|
||||
__attribute__((unused)) CONST SatFunc& satFunc = \
|
||||
*static_cast<SatFunc*>(satFunc_.get()); \
|
||||
codeToCall; \
|
||||
break; \
|
||||
} \
|
||||
};
|
||||
|
||||
SatFuncMultiplexer()
|
||||
{}
|
||||
|
||||
~SatFuncMultiplexer()
|
||||
{}
|
||||
|
||||
// this constructor not really a copy constructor, but it is
|
||||
// required to make std::unique_ptr happy
|
||||
SatFuncMultiplexer(const SatFuncMultiplexer&)
|
||||
{}
|
||||
|
||||
// this operator does not do anything and is thus not a copy operator, but it is
|
||||
// required to make std::unique_ptr happy on old compilers
|
||||
SatFuncMultiplexer& operator=(const SatFuncMultiplexer&)
|
||||
{ return *this; }
|
||||
|
||||
/*!
|
||||
* \brief Pick the correct saturation function type and initialize the object using
|
||||
* an ECL deck.
|
||||
*/
|
||||
void initFromDeck(Opm::EclipseStateConstPtr eclState,
|
||||
size_t tableIdx,
|
||||
SatFuncType satFuncType)
|
||||
{
|
||||
auto phaseUsage = phaseUsageFromDeck(eclState);
|
||||
|
||||
setSatFuncType(satFuncType);
|
||||
|
||||
OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.init(eclState,
|
||||
tableIdx,
|
||||
phaseUsage,
|
||||
/*samples=*/0), OPM_MULTIPLEXER_NON_CONST);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the type of the actually selected saturation function.
|
||||
*/
|
||||
void setSatFuncType(SatFuncType newSatFuncType)
|
||||
{
|
||||
// set the type of the saturation function
|
||||
satFuncType_ = newSatFuncType;
|
||||
|
||||
// call the default constructor for the selected saturation function
|
||||
switch (satFuncType_) {
|
||||
case Gwseg: {
|
||||
typedef SatFuncGwseg<TableType> SatFunc;
|
||||
satFunc_.reset(new SatFunc);
|
||||
break;
|
||||
}
|
||||
case Stone2: {
|
||||
typedef SatFuncStone2<TableType> SatFunc;
|
||||
satFunc_.reset(new SatFunc);
|
||||
break;
|
||||
}
|
||||
case Simple: {
|
||||
typedef SatFuncSimple<TableType> SatFunc;
|
||||
satFunc_.reset(new SatFunc);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns which saturation function is actually selected.
|
||||
*/
|
||||
SatFuncType satFuncType() const
|
||||
{ return satFuncType_; }
|
||||
|
||||
/*!
|
||||
* \brief Returns a pointer to the base saturation function.
|
||||
*
|
||||
* All underlying saturation functions derive from this class. It stores the stuff
|
||||
* which is common for all saturation functions (e.g. the endpoint scaling
|
||||
* parameters).
|
||||
*/
|
||||
const SatFuncBase<TableType>& getSatFuncBase() const
|
||||
{ return *satFunc_; }
|
||||
|
||||
SatFuncBase<TableType>& getSatFuncBase()
|
||||
{ return *satFunc_; }
|
||||
|
||||
/*!
|
||||
* \brief Return the raw "Gwseg" saturation function object.
|
||||
*
|
||||
* If the saturation function type is not "Gwseg", an assertation is triggered for
|
||||
* debug builds.
|
||||
*/
|
||||
SatFuncGwseg<TableType>& getGwseg()
|
||||
{
|
||||
assert(satFuncType_ == Gwseg);
|
||||
return *static_cast<SatFuncGwseg<TableType>*>(satFunc_.get());
|
||||
}
|
||||
|
||||
const SatFuncGwseg<TableType>& getGwseg() const
|
||||
{
|
||||
assert(satFuncType_ == Gwseg);
|
||||
return *static_cast<const SatFuncGwseg<TableType>*>(satFunc_.get());
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the raw "Simple" saturation function object.
|
||||
*
|
||||
* If the saturation function type is not "Simple", an assertation is triggered for
|
||||
* debug builds.
|
||||
*/
|
||||
SatFuncSimple<TableType>& getSimple()
|
||||
{
|
||||
assert(satFuncType_ == Simple);
|
||||
return *static_cast<SatFuncSimple<TableType>*>(satFunc_.get());
|
||||
}
|
||||
|
||||
const SatFuncSimple<TableType>& getSimple() const
|
||||
{
|
||||
assert(satFuncType_ == Simple);
|
||||
return *static_cast<const SatFuncSimple<TableType>*>(satFunc_.get());
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the raw "Stone2" saturation function object.
|
||||
*
|
||||
* If the saturation function type is not "Stone2", an assertation is triggered for
|
||||
* debug builds.
|
||||
*/
|
||||
SatFuncStone2<TableType>& getStone2()
|
||||
{
|
||||
assert(satFuncType_ == Stone2);
|
||||
return *static_cast<SatFuncStone2<TableType>*>(satFunc_.get());
|
||||
}
|
||||
|
||||
const SatFuncStone2<TableType>& getStone2() const
|
||||
{
|
||||
assert(satFuncType_ == Stone2);
|
||||
return *static_cast<const SatFuncStone2<TableType>*>(satFunc_.get());
|
||||
}
|
||||
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& fluidState, double* kr) const
|
||||
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKr(fluidState, kr), OPM_MULTIPLEXER_CONST); }
|
||||
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
|
||||
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKrDeriv(fluidState, kr, dkrds), OPM_MULTIPLEXER_CONST); }
|
||||
|
||||
template <class FluidState>
|
||||
void evalPc(const FluidState& fluidState, double* pc) const
|
||||
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPc(fluidState, pc), OPM_MULTIPLEXER_CONST); }
|
||||
|
||||
template <class FluidState>
|
||||
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
|
||||
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPcDeriv(fluidState, pc, dpcds), OPM_MULTIPLEXER_CONST); }
|
||||
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst) const
|
||||
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKr(fluidState, kr, epst), OPM_MULTIPLEXER_CONST); }
|
||||
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& fluidState, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
|
||||
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKr(fluidState, kr, epst, epst_hyst, sat_hyst), OPM_MULTIPLEXER_CONST); }
|
||||
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const
|
||||
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKrDeriv(fluidState, kr, dkrds, epst), OPM_MULTIPLEXER_CONST); }
|
||||
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
|
||||
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalKrDeriv(fluidState, kr, dkrds, epst, epst_hyst, sat_hyst), OPM_MULTIPLEXER_CONST); }
|
||||
|
||||
template <class FluidState>
|
||||
void evalPc(const FluidState& fluidState, double* pc, const EPSTransforms* epst) const
|
||||
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPc(fluidState, pc, epst), OPM_MULTIPLEXER_CONST); }
|
||||
|
||||
template <class FluidState>
|
||||
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds, const EPSTransforms* epst) const
|
||||
{ OPM_MULTIPLEXER_SATFUNC_CALL(satFunc.evalPcDeriv(fluidState, pc, dpcds, epst), OPM_MULTIPLEXER_CONST); }
|
||||
|
||||
// undefine the helper macro here because we don't need it anymore
|
||||
#undef OPM_MULTIPLEXER_SATFUNC_CALL
|
||||
#undef OPM_MULTIPLEXER_NON_CONST
|
||||
#undef OPM_MULTIPLEXER_CONST
|
||||
|
||||
private:
|
||||
SatFuncType satFuncType_;
|
||||
std::unique_ptr<SatFuncBase<TableType> > satFunc_;
|
||||
};
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_SAT_FUNC_MULTIPLEXER_HPP
|
||||
@@ -1,223 +0,0 @@
|
||||
/*
|
||||
Copyright 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/props/satfunc/SatFuncSimple.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/props/phaseUsageFromDeck.hpp>
|
||||
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <iostream>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
// SatFuncSimple.cpp can now be removed and the code below moved to SatFuncBase.cpp ...
|
||||
|
||||
template<>
|
||||
void SatFuncBase<NonuniformTableLinear<double> >::initializeTableType(NonuniformTableLinear<double> & table,
|
||||
const std::vector<double>& arg,
|
||||
const std::vector<double>& value,
|
||||
const int /* samples */)
|
||||
{
|
||||
table = NonuniformTableLinear<double>(arg, value);
|
||||
}
|
||||
|
||||
template<>
|
||||
void SatFuncBase<UniformTableLinear<double> >::initializeTableType(UniformTableLinear<double> & table,
|
||||
const std::vector<double>& arg,
|
||||
const std::vector<double>& value,
|
||||
const int samples)
|
||||
{
|
||||
buildUniformMonotoneTable(arg, value, samples, table);
|
||||
}
|
||||
|
||||
double EPSTransforms::Transform::scaleSat(double s, double s_r, double s_cr, double s_max) const
|
||||
{
|
||||
if (doNotScale) {
|
||||
return s;
|
||||
} else if (!do_3pt) { // 2-pt
|
||||
if (s <= scr) {
|
||||
return s_cr;
|
||||
} else {
|
||||
return (s >= smax) ? s_max : s_cr + (s-scr)*slope1;
|
||||
}
|
||||
} else if (s <= sr) {
|
||||
return (s <= scr) ? s_cr : s_cr+(s-scr)*slope1;
|
||||
} else {
|
||||
return (s >= smax) ? s_max : s_r+(s-sr)*slope2;
|
||||
}
|
||||
}
|
||||
|
||||
double EPSTransforms::Transform::scaleSatInv(double ss, double s_r, double s_cr, double s_max) const
|
||||
{
|
||||
if (doNotScale) {
|
||||
return ss;
|
||||
} else if (!do_3pt) { // 2-pt
|
||||
if (ss <= s_cr) {
|
||||
return scr;
|
||||
} else {
|
||||
return (ss >= s_max) ? smax : scr + (ss-s_cr)/slope1;
|
||||
}
|
||||
} else if (ss <= s_r) {
|
||||
return (ss <= s_cr) ? scr : scr+(ss-s_cr)/slope1;
|
||||
} else {
|
||||
return (ss >= s_max) ? smax : sr+(ss-s_r)/slope2;
|
||||
}
|
||||
}
|
||||
|
||||
double EPSTransforms::Transform::scaleSatDeriv(double s, double /* s_r */, double /* s_cr */, double /* s_max */) const
|
||||
{
|
||||
if (doNotScale) {
|
||||
return 1.0;
|
||||
} else if (!do_3pt) { // 2-pt
|
||||
if (s <= scr) {
|
||||
return 0.0;
|
||||
} else {
|
||||
return (s >= smax) ? 0.0 : slope1;
|
||||
}
|
||||
} else if (s <= sr) {
|
||||
return (s <= scr) ? 0.0 : slope1;
|
||||
} else {
|
||||
return (s >= smax) ? 0.0 : slope2;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double EPSTransforms::Transform::scaleSatPc(double s, double s_min, double s_max) const
|
||||
{
|
||||
if (doNotScale) {
|
||||
return s;
|
||||
} else if (s<=smin) {
|
||||
return s_min;
|
||||
} else if (s <= smax) {
|
||||
return s_min + (s-smin)*(s_max-s_min)/(smax-smin);
|
||||
} else {
|
||||
return s_max;
|
||||
}
|
||||
}
|
||||
double EPSTransforms::Transform::scaleSatDerivPc(double s, double s_min, double s_max) const
|
||||
{
|
||||
if (doNotScale) {
|
||||
return 1.0;
|
||||
} else if (s<smin) {
|
||||
return 0.0;
|
||||
} else if (s <= smax) {
|
||||
return (s_max-s_min)/(smax-smin);
|
||||
} else {
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
double EPSTransforms::Transform::scaleKr(double s, double kr, double krsr_tab) const
|
||||
{
|
||||
if (doKrCrit) {
|
||||
if (s <= scr) {
|
||||
return 0.0;
|
||||
} else if (s <= sr) {
|
||||
return kr*krSlopeCrit;
|
||||
} else if (s <= smax) {
|
||||
if (doSatInterp)
|
||||
return krsr + (s-sr)*krSlopeMax; // Note: Scaling independent of kr-value ...
|
||||
else
|
||||
return krsr + (kr-krsr_tab)*krSlopeMax;
|
||||
} else {
|
||||
return krmax;
|
||||
}
|
||||
} else if (doKrMax) {
|
||||
if (s <= scr) {
|
||||
return 0.0;
|
||||
} else if (s <= smax) {
|
||||
return kr*krSlopeMax;
|
||||
} else {
|
||||
return krmax;
|
||||
}
|
||||
} else {
|
||||
return kr;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double EPSTransforms::Transform::scaleKrDeriv(double s, double krDeriv) const
|
||||
{
|
||||
if (doKrCrit) {
|
||||
if (s < scr) {
|
||||
return 0.0;
|
||||
} else if (s <= sr) {
|
||||
return krDeriv*krSlopeCrit;
|
||||
} else if (s <= smax) {
|
||||
if (doSatInterp)
|
||||
return krSlopeMax; // Note: Scaling independent of kr-value ...
|
||||
else
|
||||
return krDeriv*krSlopeMax;
|
||||
} else {
|
||||
return 0.0;
|
||||
}
|
||||
} else if (doKrMax) {
|
||||
if (s < scr) {
|
||||
return 0.0;
|
||||
} else if (s <= smax) {
|
||||
return krDeriv*krSlopeMax;
|
||||
} else {
|
||||
return 0.0;
|
||||
}
|
||||
} else {
|
||||
if (s < scr) {
|
||||
return 0.0;
|
||||
} else if (s <= smax) {
|
||||
return krDeriv;
|
||||
} else {
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void EPSTransforms::Transform::printMe(std::ostream & out)
|
||||
{
|
||||
|
||||
out << "doNotScale: " << doNotScale << std::endl;
|
||||
out << "do_3pt: " << do_3pt << std::endl;
|
||||
out << "smin: " << smin << std::endl;
|
||||
out << "scr: " << scr << std::endl;
|
||||
out << "sr: " << sr << std::endl;
|
||||
out << "smax: " << smax << std::endl;
|
||||
out << "slope1: " << slope1 << std::endl;
|
||||
out << "slope2: " << slope2 << std::endl;
|
||||
out << "doKrMax: " << doKrMax << std::endl;
|
||||
out << "doKrCrit: " << doKrCrit << std::endl;
|
||||
out << "doSatInterp: " << doSatInterp << std::endl;
|
||||
out << "krsr: " << krsr << std::endl;
|
||||
out << "krmax: " << krmax << std::endl;
|
||||
out << "krSlopeMax: " << krSlopeMax << std::endl;
|
||||
out << "krSlopeCrit: " << krSlopeCrit << std::endl;
|
||||
}
|
||||
|
||||
void SatHyst::printMe(std::ostream & out)
|
||||
{
|
||||
out << "sg_hyst: " << sg_hyst << std::endl;
|
||||
out << "sg_shift: " << sg_shift << std::endl;
|
||||
out << "sow_hyst: " << sow_hyst << std::endl;
|
||||
out << "sow_shift: " << sow_shift << std::endl;
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
@@ -1,297 +0,0 @@
|
||||
/*
|
||||
Copyright 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#ifndef SATFUNCSIMPLE_HPP
|
||||
#define SATFUNCSIMPLE_HPP
|
||||
|
||||
#include <opm/core/props/satfunc/SatFuncBase.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
template<class TableType>
|
||||
class SatFuncSimple : public SatFuncBase<TableType>
|
||||
{
|
||||
public:
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& fluidState, double* kr) const;
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const;
|
||||
template <class FluidState>
|
||||
void evalPc(const FluidState& fluidState, double* pc) const;
|
||||
template <class FluidState>
|
||||
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const;
|
||||
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& /* fluidState */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */) const;
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& /* fluidState */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
|
||||
template <class FluidState>
|
||||
void evalPc(const FluidState& /* fluidState */, double* /* pc */, const EPSTransforms* /* epst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
|
||||
template <class FluidState>
|
||||
void evalPcDeriv(const FluidState& /* fluidState */, double* /* pc */, double* /* dpcds */, const EPSTransforms* /* epst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
|
||||
typedef SatFuncSimple<UniformTableLinear<double> > SatFuncSimpleUniform;
|
||||
typedef SatFuncSimple<NonuniformTableLinear<double> > SatFuncSimpleNonuniform;
|
||||
|
||||
template<class TableType>
|
||||
template<class FluidState>
|
||||
void SatFuncSimple<TableType>::evalKr(const FluidState& fluidState, double* kr) const
|
||||
{
|
||||
if (this->phase_usage.num_phases == 3) {
|
||||
// A simplified relative permeability model.
|
||||
double sw = fluidState.saturation(BlackoilPhases::Aqua);
|
||||
double sg = fluidState.saturation(BlackoilPhases::Vapour);
|
||||
double krw = this->krw_(sw);
|
||||
double krg = this->krg_(sg);
|
||||
double krow = this->krow_(sw + sg); // = 1 - so
|
||||
// double krog = krog_(sg); // = 1 - so - sw
|
||||
// double krocw = krocw_;
|
||||
kr[BlackoilPhases::Aqua] = krw;
|
||||
kr[BlackoilPhases::Vapour] = krg;
|
||||
kr[BlackoilPhases::Liquid] = krow;
|
||||
if (kr[BlackoilPhases::Liquid] < 0.0) {
|
||||
kr[BlackoilPhases::Liquid] = 0.0;
|
||||
}
|
||||
return;
|
||||
}
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double so = fluidState.saturation(opos);
|
||||
double krow = this->krow_(1.0-so);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double krog = this->krog_(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
}
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template<class FluidState>
|
||||
void SatFuncSimple<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
|
||||
{
|
||||
const int np = this->phase_usage.num_phases;
|
||||
std::fill(dkrds, dkrds + np*np, 0.0);
|
||||
|
||||
if (np == 3) {
|
||||
// A simplified relative permeability model.
|
||||
double sw = fluidState.saturation(BlackoilPhases::Aqua);
|
||||
double sg = fluidState.saturation(BlackoilPhases::Vapour);
|
||||
double krw = this->krw_(sw);
|
||||
double dkrww = this->krw_.derivative(sw);
|
||||
double krg = this->krg_(sg);
|
||||
double dkrgg = this->krg_.derivative(sg);
|
||||
double krow = this->krow_(sw + sg);
|
||||
double dkrow = this->krow_.derivative(sw + sg);
|
||||
// double krog = krog_(sg);
|
||||
// double dkrog = krog_.derivative(sg);
|
||||
// double krocw = krocw_;
|
||||
kr[BlackoilPhases::Aqua] = krw;
|
||||
kr[BlackoilPhases::Vapour] = krg;
|
||||
kr[BlackoilPhases::Liquid] = krow;
|
||||
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
|
||||
if (kr[BlackoilPhases::Liquid] < 0.0) {
|
||||
kr[BlackoilPhases::Liquid] = 0.0;
|
||||
}
|
||||
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
|
||||
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
|
||||
//dkrds[Liquid + Aqua*np] = dkrow;
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Liquid*np] = -dkrow;
|
||||
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = 0.0;
|
||||
//krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
|
||||
//+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
|
||||
return;
|
||||
}
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double dkrww = this->krw_.derivative(sw);
|
||||
double so = fluidState.saturation(opos);
|
||||
double krow = this->krow_(1.0-so);
|
||||
double dkrow = this->krow_.derivative(1.0-so);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
dkrds[wpos + wpos*np] = dkrww;
|
||||
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double dkrgg = this->krg_.derivative(sg);
|
||||
double krog = this->krog_(sg);
|
||||
double dkrog = this->krog_.derivative(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
dkrds[gpos + gpos*np] = dkrgg;
|
||||
dkrds[opos + gpos*np] = dkrog;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template<class FluidState>
|
||||
void SatFuncSimple<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds, const EPSTransforms* epst) const
|
||||
{
|
||||
const int np = this->phase_usage.num_phases;
|
||||
std::fill(dkrds, dkrds + np*np, 0.0);
|
||||
|
||||
if (np == 3) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
// A simplified relative permeability model.
|
||||
// Define KR(s) = scaleKr(kr(scalSat(s)))
|
||||
// Thus KR'(s) = scaleKr'(kr(scaleSat(s)))*kr'((scaleSat(s))*scaleSat'(s)
|
||||
double _sw = epst->wat.scaleSat(fluidState.saturation(wpos), 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
|
||||
double _dsdsw = epst->wat.scaleSatDeriv(fluidState.saturation(wpos), 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
|
||||
double _sg = epst->gas.scaleSat(fluidState.saturation(gpos), 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
double _dsdsg = epst->gas.scaleSatDeriv(fluidState.saturation(gpos), 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
|
||||
double _sow = epst->watoil.scaleSat(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
double _dsdsow = epst->watoil.scaleSatDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
//double _sog = epst->gasoil.scaleSat(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
//double _dsdsog = epst->gasoil.scaleSatDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
|
||||
|
||||
double krw = epst->wat.scaleKr(fluidState.saturation(wpos), this->krw_(_sw), this->krwr_);
|
||||
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(fluidState.saturation(wpos), this->krw_.derivative(_sw));
|
||||
double krg = epst->gas.scaleKr(fluidState.saturation(gpos), this->krg_(_sg), this->krgr_);
|
||||
double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(fluidState.saturation(gpos), this->krg_.derivative(_sg));
|
||||
// TODO Check the arguments to the krow- and krog-tables below...
|
||||
double krow = epst->watoil.scaleKr(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krow_(1.0-_sow-this->smin_[gpos]), this->krorw_); // ????
|
||||
double dkrow = _dsdsow*epst->watoil.scaleKrDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krow_.derivative(1.0-_sow-this->smin_[gpos])); // ????
|
||||
//double krog = epst->gasoil.scaleKr(this->krog_(1.0-_sog-this->smin_[wpos]), 1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krorg_); // ????
|
||||
//double dkrog = _dsdsog*epst->gasoil.scaleKrDeriv(1.0-fluidState.saturation(wpos)-fluidState.saturation(gpos), this->krog_.derivative(1.0-_sog-this->smin_[wpos])); // ????
|
||||
// double krocw = krocw_;
|
||||
kr[wpos] = krw;
|
||||
kr[gpos] = krg;
|
||||
kr[BlackoilPhases::Liquid] = krow;
|
||||
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
|
||||
if (kr[BlackoilPhases::Liquid] < 0.0) {
|
||||
kr[BlackoilPhases::Liquid] = 0.0;
|
||||
}
|
||||
dkrds[wpos + wpos*np] = dkrww;
|
||||
dkrds[gpos + gpos*np] = dkrgg;
|
||||
//dkrds[Liquid + Aqua*np] = dkrow;
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Liquid*np] = -dkrow;
|
||||
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
|
||||
dkrds[BlackoilPhases::Liquid + gpos*np] = 0.0;
|
||||
//krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
|
||||
//+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
|
||||
return;
|
||||
}
|
||||
OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double dkrww = this->krw_.derivative(sw);
|
||||
double so = fluidState.saturation(opos);
|
||||
double krow = this->krow_(1.0-so);
|
||||
double dkrow = this->krow_.derivative(1.0-so);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
dkrds[wpos + wpos*np] = dkrww;
|
||||
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double dkrgg = this->krg_.derivative(sg);
|
||||
double krog = this->krog_(sg);
|
||||
double dkrog = this->krog_.derivative(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
dkrds[gpos + gpos*np] = dkrgg;
|
||||
dkrds[opos + gpos*np] = dkrog;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template<class FluidState>
|
||||
void SatFuncSimple<TableType>::evalPc(const FluidState& fluidState, double* pc) const
|
||||
{
|
||||
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
pc[pos] = this->pcow_(fluidState.saturation(pos));
|
||||
}
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
pc[pos] = this->pcog_(fluidState.saturation(pos));
|
||||
}
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template<class FluidState>
|
||||
void SatFuncSimple<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
|
||||
{
|
||||
// The problem of determining three-phase capillary pressures
|
||||
// is very hard experimentally, usually one extends two-phase
|
||||
// data (as for relative permeability).
|
||||
// In our approach the derivative matrix is quite sparse, only
|
||||
// the diagonal elements corresponding to non-oil phases are
|
||||
// (potentially) nonzero.
|
||||
const int np = this->phase_usage.num_phases;
|
||||
std::fill(dpcds, dpcds + np*np, 0.0);
|
||||
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
pc[pos] = this->pcow_(fluidState.saturation(pos));
|
||||
dpcds[np*pos + pos] = this->pcow_.derivative(fluidState.saturation(pos));
|
||||
}
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
pc[pos] = this->pcog_(fluidState.saturation(pos));
|
||||
dpcds[np*pos + pos] = this->pcog_.derivative(fluidState.saturation(pos));
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
#endif // SATFUNCSIMPLE_HPP
|
||||
@@ -1,211 +0,0 @@
|
||||
/*
|
||||
Copyright 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#ifndef SATFUNCSTONE2_HPP
|
||||
#define SATFUNCSTONE2_HPP
|
||||
|
||||
#include <opm/core/props/satfunc/SatFuncBase.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
template<class TableType>
|
||||
class SatFuncStone2 : public SatFuncBase<TableType>
|
||||
{
|
||||
public:
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& fluidState, double* kr) const;
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const;
|
||||
template <class FluidState>
|
||||
void evalPc(const FluidState& fluidState, double* pc) const;
|
||||
template <class FluidState>
|
||||
void evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const;
|
||||
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
|
||||
template <class FluidState>
|
||||
void evalKr(const FluidState& /* fluidState */, double* /* kr */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& /* fluidState */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
|
||||
template <class FluidState>
|
||||
void evalKrDeriv(const FluidState& /* fluidState */, double* /* kr */, double* /* dkrds */, const EPSTransforms* /* epst */, const EPSTransforms* /* epst_hyst */, const SatHyst* /* sat_hyst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
|
||||
template <class FluidState>
|
||||
void evalPc(const FluidState& /* fluidState */, double* /* pc */, const EPSTransforms* /* epst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
|
||||
template <class FluidState>
|
||||
void evalPcDeriv(const FluidState& /* fluidState */, double* /* pc */, double* /* dpcds */, const EPSTransforms* /* epst */) const
|
||||
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
|
||||
typedef SatFuncStone2<UniformTableLinear<double> > SatFuncStone2Uniform;
|
||||
typedef SatFuncStone2<NonuniformTableLinear<double> > SatFuncStone2Nonuniform;
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncStone2<TableType>::evalKr(const FluidState& fluidState, double* kr) const
|
||||
{
|
||||
if (this->phase_usage.num_phases == 3) {
|
||||
// Stone-II relative permeability model.
|
||||
double sw = fluidState.saturation(BlackoilPhases::Aqua);
|
||||
double sg = fluidState.saturation(BlackoilPhases::Vapour);
|
||||
double krw = this->krw_(sw);
|
||||
double krg = this->krg_(sg);
|
||||
double krow = this->krow_(sw + sg); // = 1 - so
|
||||
double krog = this->krog_(sg); // = 1 - so - sw
|
||||
double krocw = this->krocw_;
|
||||
kr[BlackoilPhases::Aqua] = krw;
|
||||
kr[BlackoilPhases::Vapour] = krg;
|
||||
kr[BlackoilPhases::Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
|
||||
if (kr[BlackoilPhases::Liquid] < 0.0) {
|
||||
kr[BlackoilPhases::Liquid] = 0.0;
|
||||
}
|
||||
return;
|
||||
}
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double krow = this->krow_(sw);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double krog = this->krog_(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
}
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncStone2<TableType>::evalKrDeriv(const FluidState& fluidState, double* kr, double* dkrds) const
|
||||
{
|
||||
const int np = this->phase_usage.num_phases;
|
||||
std::fill(dkrds, dkrds + np*np, 0.0);
|
||||
|
||||
if (np == 3) {
|
||||
// Stone-II relative permeability model.
|
||||
double sw = fluidState.saturation(BlackoilPhases::Aqua);
|
||||
double sg = fluidState.saturation(BlackoilPhases::Vapour);
|
||||
double krw = this->krw_(sw);
|
||||
double dkrww = this->krw_.derivative(sw);
|
||||
double krg = this->krg_(sg);
|
||||
double dkrgg = this->krg_.derivative(sg);
|
||||
double krow = this->krow_(sw + sg);
|
||||
double dkrow = this->krow_.derivative(sw + sg);
|
||||
double krog = this->krog_(sg);
|
||||
double dkrog = this->krog_.derivative(sg);
|
||||
double krocw = this->krocw_;
|
||||
kr[BlackoilPhases::Aqua] = krw;
|
||||
kr[BlackoilPhases::Vapour] = krg;
|
||||
kr[BlackoilPhases::Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
|
||||
if (kr[BlackoilPhases::Liquid] < 0.0) {
|
||||
kr[BlackoilPhases::Liquid] = 0.0;
|
||||
}
|
||||
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
|
||||
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
|
||||
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
|
||||
+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
|
||||
return;
|
||||
}
|
||||
// We have a two-phase situation. We know that oil is active.
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sw = fluidState.saturation(wpos);
|
||||
double krw = this->krw_(sw);
|
||||
double dkrww = this->krw_.derivative(sw);
|
||||
double krow = this->krow_(sw);
|
||||
double dkrow = this->krow_.derivative(sw);
|
||||
kr[wpos] = krw;
|
||||
kr[opos] = krow;
|
||||
dkrds[wpos + wpos*np] = dkrww;
|
||||
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
|
||||
} else {
|
||||
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
|
||||
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
|
||||
double sg = fluidState.saturation(gpos);
|
||||
double krg = this->krg_(sg);
|
||||
double dkrgg = this->krg_.derivative(sg);
|
||||
double krog = this->krog_(sg);
|
||||
double dkrog = this->krog_.derivative(sg);
|
||||
kr[gpos] = krg;
|
||||
kr[opos] = krog;
|
||||
dkrds[gpos + gpos*np] = dkrgg;
|
||||
dkrds[opos + gpos*np] = dkrog;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncStone2<TableType>::evalPc(const FluidState& fluidState, double* pc) const
|
||||
{
|
||||
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
pc[pos] = this->pcow_(fluidState.saturation(pos));
|
||||
}
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
pc[pos] = this->pcog_(fluidState.saturation(pos));
|
||||
}
|
||||
}
|
||||
|
||||
template<class TableType>
|
||||
template <class FluidState>
|
||||
void SatFuncStone2<TableType>::evalPcDeriv(const FluidState& fluidState, double* pc, double* dpcds) const
|
||||
{
|
||||
// The problem of determining three-phase capillary pressures
|
||||
// is very hard experimentally, usually one extends two-phase
|
||||
// data (as for relative permeability).
|
||||
// In our approach the derivative matrix is quite sparse, only
|
||||
// the diagonal elements corresponding to non-oil phases are
|
||||
// (potentially) nonzero.
|
||||
const int np = this->phase_usage.num_phases;
|
||||
std::fill(dpcds, dpcds + np*np, 0.0);
|
||||
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
|
||||
pc[pos] = this->pcow_(fluidState.saturation(pos));
|
||||
dpcds[np*pos + pos] = this->pcow_.derivative(fluidState.saturation(pos));
|
||||
}
|
||||
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
|
||||
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
|
||||
pc[pos] = this->pcog_(fluidState.saturation(pos));
|
||||
dpcds[np*pos + pos] = this->pcog_.derivative(fluidState.saturation(pos));
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
#endif // SATFUNCSTONE2_HPP
|
||||
269
opm/core/props/satfunc/SaturationPropsFromDeck.cpp
Normal file
269
opm/core/props/satfunc/SaturationPropsFromDeck.cpp
Normal file
@@ -0,0 +1,269 @@
|
||||
/*
|
||||
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/utility/UniformTableLinear.hpp>
|
||||
#include <opm/core/utility/NonuniformTableLinear.hpp>
|
||||
#include <opm/core/grid/GridHelpers.hpp>
|
||||
#include <opm/core/simulator/ExplicitArraysFluidState.hpp>
|
||||
#include <opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/EndscaleWrapper.hpp>
|
||||
#include <opm/parser/eclipse/Utility/ScalecrsWrapper.hpp>
|
||||
|
||||
#include <iostream>
|
||||
#include <map>
|
||||
|
||||
#include "SaturationPropsFromDeck.hpp"
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
|
||||
// ----------- Methods of SaturationPropsFromDeck ---------
|
||||
|
||||
|
||||
/// Default constructor.
|
||||
SaturationPropsFromDeck::SaturationPropsFromDeck()
|
||||
{
|
||||
}
|
||||
|
||||
/// Initialize from deck.
|
||||
void SaturationPropsFromDeck::init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const UnstructuredGrid& grid)
|
||||
{
|
||||
this->init(deck, eclipseState, materialLawManager, grid.number_of_cells,
|
||||
grid.global_cell, grid.cell_centroids,
|
||||
grid.dimensions);
|
||||
}
|
||||
|
||||
/// Initialize from deck.
|
||||
void SaturationPropsFromDeck::init(const PhaseUsage &phaseUsage,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager)
|
||||
{
|
||||
phaseUsage_ = phaseUsage;
|
||||
materialLawManager_ = materialLawManager;
|
||||
}
|
||||
|
||||
/// \return P, the number of phases.
|
||||
int SaturationPropsFromDeck::numPhases() const
|
||||
{
|
||||
return phaseUsage_.num_phases;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/// Relative permeability.
|
||||
/// \param[in] n Number of data points.
|
||||
/// \param[in] s Array of nP saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the s values.
|
||||
/// \param[out] kr Array of nP relperm values, array must be valid before calling.
|
||||
/// \param[out] dkrds If non-null: array of nP^2 relperm derivative values,
|
||||
/// array must be valid before calling.
|
||||
/// 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 SaturationPropsFromDeck::relperm(const int n,
|
||||
const double* s,
|
||||
const int* cells,
|
||||
double* kr,
|
||||
double* dkrds) const
|
||||
{
|
||||
assert(cells != 0);
|
||||
|
||||
const int np = BlackoilPhases::MaxNumPhases;
|
||||
if (dkrds) {
|
||||
ExplicitArraysSatDerivativesFluidState fluidState(phaseUsage_);
|
||||
fluidState.setSaturationArray(s);
|
||||
|
||||
typedef ExplicitArraysSatDerivativesFluidState::Evaluation Evaluation;
|
||||
Evaluation relativePerms[BlackoilPhases::MaxNumPhases];
|
||||
for (int i = 0; i < n; ++i) {
|
||||
fluidState.setIndex(i);
|
||||
const auto& params = materialLawManager_->materialLawParams(cells[i]);
|
||||
MaterialLaw::relativePermeabilities(relativePerms, params, fluidState);
|
||||
|
||||
// copy the values calculated using opm-material to the target arrays
|
||||
for (int krPhaseIdx = 0; krPhaseIdx < np; ++krPhaseIdx) {
|
||||
kr[np*i + krPhaseIdx] = relativePerms[krPhaseIdx].value;
|
||||
|
||||
for (int satPhaseIdx = 0; satPhaseIdx < np; ++satPhaseIdx)
|
||||
dkrds[np*np*i + satPhaseIdx*np + krPhaseIdx] = relativePerms[krPhaseIdx].derivatives[satPhaseIdx];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
ExplicitArraysFluidState fluidState(phaseUsage_);
|
||||
fluidState.setSaturationArray(s);
|
||||
|
||||
double relativePerms[BlackoilPhases::MaxNumPhases];
|
||||
for (int i = 0; i < n; ++i) {
|
||||
fluidState.setIndex(i);
|
||||
const auto& params = materialLawManager_->materialLawParams(cells[i]);
|
||||
MaterialLaw::relativePermeabilities(relativePerms, params, fluidState);
|
||||
|
||||
// copy the values calculated using opm-material to the target arrays
|
||||
for (int krPhaseIdx = 0; krPhaseIdx < np; ++krPhaseIdx) {
|
||||
kr[np*i + krPhaseIdx] = relativePerms[krPhaseIdx];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/// Capillary pressure.
|
||||
/// \param[in] n Number of data points.
|
||||
/// \param[in] s Array of nP saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the s values.
|
||||
/// \param[out] pc Array of nP capillary pressure values, array must be valid before calling.
|
||||
/// \param[out] dpcds If non-null: array of nP^2 derivative values,
|
||||
/// array must be valid before calling.
|
||||
/// 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 SaturationPropsFromDeck::capPress(const int n,
|
||||
const double* s,
|
||||
const int* cells,
|
||||
double* pc,
|
||||
double* dpcds) const
|
||||
{
|
||||
assert(cells != 0);
|
||||
|
||||
const int np = BlackoilPhases::MaxNumPhases;
|
||||
if (dpcds) {
|
||||
ExplicitArraysSatDerivativesFluidState fluidState(phaseUsage_);
|
||||
typedef ExplicitArraysSatDerivativesFluidState::Evaluation Evaluation;
|
||||
fluidState.setSaturationArray(s);
|
||||
|
||||
Evaluation capillaryPressures[BlackoilPhases::MaxNumPhases];
|
||||
for (int i = 0; i < n; ++i) {
|
||||
fluidState.setIndex(i);
|
||||
const auto& params = materialLawManager_->materialLawParams(cells[i]);
|
||||
MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState);
|
||||
|
||||
// copy the values calculated using opm-material to the target arrays
|
||||
for (int pcPhaseIdx = 0; pcPhaseIdx < np; ++pcPhaseIdx) {
|
||||
double sign = (pcPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
|
||||
pc[np*i + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx].value;
|
||||
|
||||
for (int satPhaseIdx = 0; satPhaseIdx < np; ++satPhaseIdx)
|
||||
dpcds[np*np*i + satPhaseIdx*np + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx].derivatives[satPhaseIdx];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
ExplicitArraysFluidState fluidState(phaseUsage_);
|
||||
fluidState.setSaturationArray(s);
|
||||
|
||||
double capillaryPressures[BlackoilPhases::MaxNumPhases];
|
||||
for (int i = 0; i < n; ++i) {
|
||||
fluidState.setIndex(i);
|
||||
const auto& params = materialLawManager_->materialLawParams(cells[i]);
|
||||
MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState);
|
||||
|
||||
// copy the values calculated using opm-material to the target arrays
|
||||
for (int pcPhaseIdx = 0; pcPhaseIdx < np; ++pcPhaseIdx) {
|
||||
double sign = (pcPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
|
||||
pc[np*i + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// 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.
|
||||
void SaturationPropsFromDeck::satRange(const int n,
|
||||
const int* cells,
|
||||
double* smin,
|
||||
double* smax) const
|
||||
{
|
||||
int wpos = phaseUsage_.phase_pos[BlackoilPhases::Aqua];
|
||||
int gpos = phaseUsage_.phase_pos[BlackoilPhases::Vapour];
|
||||
int opos = phaseUsage_.phase_pos[BlackoilPhases::Liquid];
|
||||
|
||||
const int np = phaseUsage_.num_phases;
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const auto& scaledDrainageInfo =
|
||||
materialLawManager_->oilWaterScaledEpsInfoDrainage(cells[i]);
|
||||
|
||||
if (phaseUsage_.phase_used[BlackoilPhases::Aqua]) {
|
||||
smin[np*i + wpos] = scaledDrainageInfo.Swl;
|
||||
smax[np*i + wpos] = scaledDrainageInfo.Swu;
|
||||
}
|
||||
|
||||
if (phaseUsage_.phase_used[BlackoilPhases::Vapour]) {
|
||||
smin[np*i + gpos] = scaledDrainageInfo.Sgl;
|
||||
smax[np*i + gpos] = scaledDrainageInfo.Sgu;
|
||||
}
|
||||
|
||||
if (phaseUsage_.phase_used[BlackoilPhases::Liquid]) {
|
||||
smin[np*i + opos] = 1.0;
|
||||
smax[np*i + opos] = 1.0;
|
||||
if (phaseUsage_.phase_used[BlackoilPhases::Aqua]) {
|
||||
smin[np*i + opos] -= smax[np*i + wpos];
|
||||
smax[np*i + opos] -= smin[np*i + wpos];
|
||||
}
|
||||
if (phaseUsage_.phase_used[BlackoilPhases::Vapour]) {
|
||||
smin[np*i + opos] -= smax[np*i + gpos];
|
||||
smax[np*i + opos] -= smin[np*i + gpos];
|
||||
}
|
||||
smin[np*i + opos] = std::max(0.0, smin[np*i + opos]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Update saturation state for the hysteresis tracking
|
||||
/// \param[in] n Number of data points.
|
||||
/// \param[in] s Array of nP saturation values.
|
||||
void SaturationPropsFromDeck::updateSatHyst(const int n,
|
||||
const int* cells,
|
||||
const double* s)
|
||||
{
|
||||
assert(cells != 0);
|
||||
|
||||
if (materialLawManager_->enableHysteresis()) {
|
||||
ExplicitArraysFluidState fluidState(phaseUsage_);
|
||||
fluidState.setSaturationArray(s);
|
||||
|
||||
for (int i = 0; i < n; ++i) {
|
||||
fluidState.setIndex(i);
|
||||
materialLawManager_->updateHysteresis(fluidState, cells[i]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// 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 SaturationPropsFromDeck::swatInitScaling(const int cell,
|
||||
const double pcow,
|
||||
double& swat)
|
||||
{
|
||||
swat = materialLawManager_->applySwatinit(cell, pcow, swat);
|
||||
}
|
||||
} // namespace Opm
|
||||
@@ -23,11 +23,14 @@
|
||||
#include <opm/core/props/satfunc/SaturationPropsInterface.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/props/satfunc/SatFuncMultiplexer.hpp>
|
||||
#include <opm/core/props/phaseUsageFromDeck.hpp>
|
||||
#include <opm/core/grid.h>
|
||||
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
||||
|
||||
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
struct UnstructuredGrid;
|
||||
@@ -38,17 +41,32 @@ namespace Opm
|
||||
class SaturationPropsFromDeck : public SaturationPropsInterface
|
||||
{
|
||||
public:
|
||||
typedef Opm::ThreePhaseMaterialTraits<double,
|
||||
/*wettingPhaseIdx=*/BlackoilPhases::Aqua,
|
||||
/*nonWettingPhaseIdx=*/BlackoilPhases::Liquid,
|
||||
/*gasPhaseIdx=*/BlackoilPhases::Vapour> MaterialTraits;
|
||||
typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
|
||||
typedef MaterialLawManager::MaterialLaw MaterialLaw;
|
||||
typedef MaterialLawManager::MaterialLawParams MaterialLawParams;
|
||||
|
||||
/// Default constructor.
|
||||
inline SaturationPropsFromDeck();
|
||||
SaturationPropsFromDeck();
|
||||
|
||||
/// Initialize from a MaterialLawManager object and a compressed to cartesian cell index map.
|
||||
/// \param[in] materialLawManager An initialized MaterialLawManager object
|
||||
void init(const PhaseUsage &phaseUsage,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager);
|
||||
|
||||
|
||||
/// 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.
|
||||
inline void init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
const UnstructuredGrid& grid);
|
||||
void init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const UnstructuredGrid& grid);
|
||||
|
||||
/// Initialize from deck and grid.
|
||||
/// \param[in] deck Deck input parser
|
||||
@@ -62,15 +80,19 @@ namespace Opm
|
||||
/// \param[in] begin_cell_centroids Pointer to the first cell_centroid of the grid.
|
||||
/// \param[in] dimensions The dimensions of the grid.
|
||||
template<class T>
|
||||
inline void init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const T& begin_cell_centroids,
|
||||
int dimensions);
|
||||
void init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const T& begin_cell_centroids,
|
||||
int dimensions)
|
||||
{
|
||||
init(Opm::phaseUsageFromDeck(deck), materialLawManager);
|
||||
}
|
||||
|
||||
/// \return P, the number of phases.
|
||||
inline int numPhases() const;
|
||||
int numPhases() const;
|
||||
|
||||
/// Relative permeability.
|
||||
/// \param[in] n Number of data points.
|
||||
@@ -81,11 +103,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 ...)
|
||||
inline void relperm(const int n,
|
||||
const double* s,
|
||||
const int* cells,
|
||||
double* kr,
|
||||
double* dkrds) const;
|
||||
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.
|
||||
@@ -96,108 +118,43 @@ 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 ...)
|
||||
inline void capPress(const int n,
|
||||
const double* s,
|
||||
const int* cells,
|
||||
double* pc,
|
||||
double* dpcds) const;
|
||||
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.
|
||||
inline void satRange(const int n,
|
||||
const int* cells,
|
||||
double* smin,
|
||||
double* smax) const;
|
||||
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.
|
||||
inline void updateSatHyst(const int n,
|
||||
const int* cells,
|
||||
const double* s);
|
||||
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.
|
||||
inline void swatInitScaling(const int cell,
|
||||
const double pcow,
|
||||
double & swat);
|
||||
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_;
|
||||
typedef Opm::SatFuncMultiplexer SatFunc;
|
||||
std::vector<SatFunc> satfunc_;
|
||||
std::vector<int> cell_to_func_; // = SATNUM - 1
|
||||
std::vector<int> cell_to_func_imb_;
|
||||
|
||||
bool do_eps_; // ENDSCALE is active
|
||||
bool do_3pt_; // SCALECRS: YES~true NO~false
|
||||
bool do_hyst_; // Keywords ISWL etc detected
|
||||
std::vector<EPSTransforms> eps_transf_;
|
||||
std::vector<EPSTransforms> eps_transf_hyst_;
|
||||
std::vector<SatHyst> sat_hyst_;
|
||||
|
||||
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,
|
||||
const std::vector<std::string>& eps_kw,
|
||||
std::vector<EPSTransforms>& eps_transf);
|
||||
template<class T>
|
||||
void initEPSKey(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const T& begin_cell_centroids,
|
||||
int dimensions,
|
||||
const std::string& keyword,
|
||||
std::vector<double>& scaleparam);
|
||||
void initEPSParam(const int cell,
|
||||
EPSTransforms::Transform& data,
|
||||
const bool oil,
|
||||
const double sl_tab,
|
||||
const double scr_tab,
|
||||
const double su_tab,
|
||||
const double sxcr_tab,
|
||||
const double s0_tab,
|
||||
const double krsr_tab,
|
||||
const double krmax_tab,
|
||||
const double pcmax_tab,
|
||||
const std::vector<double>& sl,
|
||||
const std::vector<double>& scr,
|
||||
const std::vector<double>& su,
|
||||
const std::vector<double>& sxcr,
|
||||
const std::vector<double>& s0,
|
||||
const std::vector<double>& krsr,
|
||||
const std::vector<double>& krmax,
|
||||
const std::vector<double>& pcmax);
|
||||
|
||||
bool columnIsMasked_(Opm::DeckConstPtr deck,
|
||||
const std::string& keywordName,
|
||||
int columnIdx)
|
||||
{ return deck->getKeyword(keywordName)->getRecord(columnIdx)->getItem(0)->getSIDouble(0) != -1.0; }
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager_;
|
||||
PhaseUsage phaseUsage_;
|
||||
};
|
||||
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
#include <opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp>
|
||||
|
||||
|
||||
#endif // OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
|
||||
|
||||
@@ -1,864 +0,0 @@
|
||||
/*
|
||||
Copyright 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED
|
||||
#define OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED
|
||||
|
||||
|
||||
#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>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/EndscaleWrapper.hpp>
|
||||
#include <opm/parser/eclipse/Utility/ScalecrsWrapper.hpp>
|
||||
|
||||
#include <iostream>
|
||||
#include <map>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
|
||||
// ----------- Methods of SaturationPropsFromDeck ---------
|
||||
|
||||
|
||||
/// Default constructor.
|
||||
inline
|
||||
SaturationPropsFromDeck::SaturationPropsFromDeck()
|
||||
{
|
||||
}
|
||||
|
||||
/// Initialize from deck.
|
||||
inline
|
||||
void SaturationPropsFromDeck::init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
const UnstructuredGrid& grid)
|
||||
{
|
||||
this->init(deck, eclipseState, grid.number_of_cells,
|
||||
grid.global_cell, grid.cell_centroids,
|
||||
grid.dimensions);
|
||||
}
|
||||
|
||||
/// Initialize from deck.
|
||||
template<class T>
|
||||
void SaturationPropsFromDeck::init(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const T& begin_cell_centroids,
|
||||
int dimensions)
|
||||
{
|
||||
phase_usage_ = phaseUsageFromDeck(deck);
|
||||
|
||||
// Extract input data.
|
||||
// Oil phase should be active.
|
||||
if (!phase_usage_.phase_used[BlackoilPhases::Liquid]) {
|
||||
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- oil phase must be active.");
|
||||
}
|
||||
|
||||
// Check SATOPTS status
|
||||
bool hysteresis_switch = false;
|
||||
if (deck->hasKeyword("SATOPTS")) {
|
||||
const std::vector<std::string>& satopts = deck->getKeyword("SATOPTS")->getStringData();
|
||||
for (size_t i = 0; i < satopts.size(); ++i) {
|
||||
if (satopts[i] == std::string("HYSTER")) {
|
||||
hysteresis_switch = true;
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Keyword SATOPTS: Switch " << satopts[i] << " not supported. ");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 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;
|
||||
const int* gc = global_cell;
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
const int deck_pos = (gc == NULL) ? cell : gc[cell];
|
||||
cell_to_func_[cell] = satnum[deck_pos] - 1;
|
||||
}
|
||||
}
|
||||
|
||||
// Find number of tables, check for consistency.
|
||||
enum { Uninitialized = -1 };
|
||||
int num_tables = Uninitialized;
|
||||
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[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);
|
||||
}
|
||||
if (num_tables == Uninitialized) {
|
||||
num_tables = num_sgof_tables;
|
||||
} else if (num_tables != num_sgof_tables) {
|
||||
OPM_THROW(std::runtime_error, "Inconsistent number of tables in SWOF and SGOF.");
|
||||
}
|
||||
}
|
||||
|
||||
// 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")) {
|
||||
const int& relative_perm_hyst = deck->getKeyword("EHYSTR")->getRecord(0)->getItem(1)->getInt(0);
|
||||
const std::string& limiting_hyst_flag = deck->getKeyword("EHYSTR")->getRecord(0)->getItem(4)->getString(0);
|
||||
if (relative_perm_hyst != int(0)) {
|
||||
OPM_THROW(std::runtime_error, "Keyword EHYSTR, item 2: Flag '" << relative_perm_hyst << "' found, only '0' is supported. ");
|
||||
}
|
||||
if (limiting_hyst_flag != std::string("KR")) {
|
||||
OPM_THROW(std::runtime_error, "Keyword EHYSTR, item 5: Flag '" << limiting_hyst_flag << "' found, only 'KR' is supported. ");
|
||||
}
|
||||
if ( ! deck->hasKeyword("ENDSCALE")) {
|
||||
// TODO When use of IMBNUM is implemented, this constraint will be lifted.
|
||||
OPM_THROW(std::runtime_error, "Currently hysteris effects is only available through endpoint scaling.");
|
||||
}
|
||||
do_hyst_ = true;
|
||||
} else if (hysteresis_switch) {
|
||||
OPM_THROW(std::runtime_error, "Switch HYSTER of keyword SATOPTS is active, but keyword EHYSTR not found.");
|
||||
} else if (deck->hasKeyword("EHYSTR")) {
|
||||
OPM_THROW(std::runtime_error, "Found keyword EHYSTR, but switch HYSTER of keyword SATOPTS is not set.");
|
||||
}
|
||||
|
||||
// Saturation table scaling
|
||||
do_eps_ = false;
|
||||
do_3pt_ = false;
|
||||
if (deck->hasKeyword("ENDSCALE")) {
|
||||
Opm::EndscaleWrapper endscale(deck->getKeyword("ENDSCALE"));
|
||||
if (endscale.directionSwitch() != std::string("NODIR")) {
|
||||
OPM_THROW(std::runtime_error,
|
||||
"SaturationPropsFromDeck::init() -- ENDSCALE: "
|
||||
"Currently only 'NODIR' accepted.");
|
||||
}
|
||||
if (!endscale.isReversible()) {
|
||||
OPM_THROW(std::runtime_error,
|
||||
"SaturationPropsFromDeck::init() -- ENDSCALE: "
|
||||
"Currently only 'REVERS' accepted.");
|
||||
}
|
||||
if (deck->hasKeyword("SCALECRS")) {
|
||||
Opm::ScalecrsWrapper scalecrs(deck->getKeyword("SCALECRS"));
|
||||
if (scalecrs.isEnabled()) {
|
||||
do_3pt_ = true;
|
||||
}
|
||||
}
|
||||
do_eps_ = true;
|
||||
|
||||
// Make a consistency check of ENDNUM: #regions = NTENDP (ENDSCALE::3, TABDIMS::8)...
|
||||
if (deck->hasKeyword("ENDNUM")) {
|
||||
const std::vector<int>& endnum = deck->getKeyword("ENDNUM")->getIntData();
|
||||
int endnum_regions = *std::max_element(endnum.begin(), endnum.end());
|
||||
if (endnum_regions > endscale.numEndscaleTables()) {
|
||||
OPM_THROW(std::runtime_error,
|
||||
"ENDNUM: Found " << endnum_regions <<
|
||||
" regions. Maximum allowed is " << endscale.numEndscaleTables() <<
|
||||
" (confer item 3 of keyword ENDSCALE)."); // TODO See also item 8 of TABDIMS ...
|
||||
}
|
||||
}
|
||||
// TODO: ENPTVD/ENKRVD: Too few tables gives a cryptical message from parser,
|
||||
// superfluous tables are ignored by the parser without any warning ...
|
||||
|
||||
const std::vector<std::string> eps_kw{"SWL", "SWU", "SWCR", "SGL", "SGU", "SGCR", "SOWCR",
|
||||
"SOGCR", "KRW", "KRG", "KRO", "KRWR", "KRGR", "KRORW", "KRORG", "PCW", "PCG"};
|
||||
eps_transf_.resize(number_of_cells);
|
||||
initEPS(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroids,
|
||||
dimensions, eps_kw, eps_transf_);
|
||||
|
||||
if (do_hyst_) {
|
||||
if (deck->hasKeyword("KRW")
|
||||
|| deck->hasKeyword("KRG")
|
||||
|| deck->hasKeyword("KRO")
|
||||
|| deck->hasKeyword("KRWR")
|
||||
|| deck->hasKeyword("KRGR")
|
||||
|| deck->hasKeyword("KRORW")
|
||||
|| deck->hasKeyword("KRORG")
|
||||
|| deck->hasKeyword("ENKRVD")
|
||||
|| deck->hasKeyword("IKRG")
|
||||
|| deck->hasKeyword("IKRO")
|
||||
|| deck->hasKeyword("IKRWR")
|
||||
|| deck->hasKeyword("IKRGR")
|
||||
|| deck->hasKeyword("IKRORW")
|
||||
|| deck->hasKeyword("IKRORG") ) {
|
||||
OPM_THROW(std::runtime_error,"Currently hysteresis and relperm value scaling cannot be combined.");
|
||||
}
|
||||
|
||||
if (deck->hasKeyword("IMBNUM")) {
|
||||
const std::vector<int>& imbnum = deck->getKeyword("IMBNUM")->getIntData();
|
||||
int imbnum_regions = *std::max_element(imbnum.begin(), imbnum.end());
|
||||
if (imbnum_regions > num_tables) {
|
||||
OPM_THROW(std::runtime_error,
|
||||
"IMBNUM: Found " << imbnum_regions <<
|
||||
" regions. Maximum allowed is " << num_tables <<
|
||||
" (number of tables provided by SWOF/SGOF).");
|
||||
}
|
||||
const int num_cells = number_of_cells;
|
||||
cell_to_func_imb_.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];
|
||||
cell_to_func_imb_[cell] = imbnum[deck_pos] - 1;
|
||||
}
|
||||
// TODO: Make actual use of IMBNUM. For now we just consider the imbibition curve
|
||||
// to be a scaled version of the drainage curve (confer Norne model).
|
||||
}
|
||||
|
||||
const std::vector<std::string> eps_i_kw{"ISWL", "ISWU", "ISWCR", "ISGL", "ISGU", "ISGCR", "ISOWCR",
|
||||
"ISOGCR", "IKRW", "IKRG", "IKRO", "IKRWR", "IKRGR", "IKRORW", "IKRORG", "IPCW", "IPCG"};
|
||||
eps_transf_hyst_.resize(number_of_cells);
|
||||
sat_hyst_.resize(number_of_cells);
|
||||
initEPS(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroids,
|
||||
dimensions, eps_i_kw, eps_transf_hyst_);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/// \return P, the number of phases.
|
||||
inline
|
||||
int SaturationPropsFromDeck::numPhases() const
|
||||
{
|
||||
return phase_usage_.num_phases;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/// Relative permeability.
|
||||
/// \param[in] n Number of data points.
|
||||
/// \param[in] s Array of nP saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the s values.
|
||||
/// \param[out] kr Array of nP relperm values, array must be valid before calling.
|
||||
/// \param[out] dkrds If non-null: array of nP^2 relperm derivative values,
|
||||
/// array must be valid before calling.
|
||||
/// 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 ...)
|
||||
inline
|
||||
void SaturationPropsFromDeck::relperm(const int n,
|
||||
const double* s,
|
||||
const int* cells,
|
||||
double* kr,
|
||||
double* dkrds) const
|
||||
{
|
||||
assert(cells != 0);
|
||||
|
||||
const int np = phase_usage_.num_phases;
|
||||
ExplicitArraysFluidState fluidState(np);
|
||||
fluidState.setSaturationArray(s);
|
||||
|
||||
if (dkrds) {
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
fluidState.setIndex(i);
|
||||
if (do_hyst_) {
|
||||
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_) {
|
||||
satfunc_[cell_to_func_[cells[i]]].evalKrDeriv(fluidState, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]]));
|
||||
} else {
|
||||
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) {
|
||||
fluidState.setIndex(i);
|
||||
if (do_hyst_) {
|
||||
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_) {
|
||||
satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i, &(eps_transf_[cells[i]]));
|
||||
} else {
|
||||
satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/// Capillary pressure.
|
||||
/// \param[in] n Number of data points.
|
||||
/// \param[in] s Array of nP saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the s values.
|
||||
/// \param[out] pc Array of nP capillary pressure values, array must be valid before calling.
|
||||
/// \param[out] dpcds If non-null: array of nP^2 derivative values,
|
||||
/// array must be valid before calling.
|
||||
/// 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 ...)
|
||||
inline
|
||||
void SaturationPropsFromDeck::capPress(const int n,
|
||||
const double* s,
|
||||
const int* cells,
|
||||
double* pc,
|
||||
double* dpcds) const
|
||||
{
|
||||
assert(cells != 0);
|
||||
|
||||
const int np = phase_usage_.num_phases;
|
||||
ExplicitArraysFluidState fluidState(np);
|
||||
fluidState.setSaturationArray(s);
|
||||
|
||||
if (dpcds) {
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
fluidState.setIndex(i);
|
||||
if (do_eps_) {
|
||||
satfunc_[cell_to_func_[cells[i]]].evalPcDeriv(fluidState, pc + np*i, dpcds + np*np*i, &(eps_transf_[cells[i]]));
|
||||
} else {
|
||||
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_) {
|
||||
satfunc_[cell_to_func_[cells[i]]].evalPc(fluidState, pc + np*i, &(eps_transf_[cells[i]]));
|
||||
} else {
|
||||
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.
|
||||
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;
|
||||
|
||||
if (do_eps_) {
|
||||
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
|
||||
const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid];
|
||||
const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
|
||||
|
||||
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 p = 0; p < np; ++p) {
|
||||
smin[np*cellIdx + p] = satFunc.smin_[p];
|
||||
smax[np*cellIdx + p] = satFunc.smax_[p];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// Update saturation state for the hysteresis tracking
|
||||
/// \param[in] n Number of data points.
|
||||
/// \param[in] s Array of nP saturation values.
|
||||
inline
|
||||
void SaturationPropsFromDeck::updateSatHyst(const int n,
|
||||
const int* cells,
|
||||
const double* s)
|
||||
{
|
||||
assert(cells != 0);
|
||||
|
||||
const int np = phase_usage_.num_phases;
|
||||
if (do_hyst_) {
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
satfunc_[cell_to_func_[cells[i]]].getSatFuncBase().updateSatHyst(s + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]]));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// 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.
|
||||
inline
|
||||
void SaturationPropsFromDeck::swatInitScaling(const int cell,
|
||||
const double pcow,
|
||||
double& swat)
|
||||
{
|
||||
if (phase_usage_.phase_used[BlackoilPhases::Aqua]) {
|
||||
const double pc_low_threshold = 1.0e-8;
|
||||
// TODO: Mixed wettability systems - see ecl kw OPTIONS switch 74
|
||||
if (swat <= eps_transf_[cell].wat.smin) {
|
||||
swat = eps_transf_[cell].wat.smin;
|
||||
} else if (pcow < pc_low_threshold) {
|
||||
swat = eps_transf_[cell].wat.smax;
|
||||
} else {
|
||||
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
|
||||
const int max_np = BlackoilPhases::MaxNumPhases;
|
||||
double s[max_np] = { 0.0 };
|
||||
s[wpos] = swat;
|
||||
ExplicitArraysFluidState fluidState(phase_usage_.num_phases);
|
||||
fluidState.setSaturationArray(s);
|
||||
fluidState.setIndex(0);
|
||||
double pc[max_np] = { 0.0 };
|
||||
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];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "swatInitScaling: no water phase! ");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Initialize saturation scaling parameters
|
||||
template<class T>
|
||||
void SaturationPropsFromDeck::initEPS(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const T& begin_cell_centroid,
|
||||
int dimensions,
|
||||
const std::vector<std::string>& eps_kw,
|
||||
std::vector<EPSTransforms>& eps_transf)
|
||||
{
|
||||
std::vector<std::vector<double> > eps_vec(eps_kw.size());
|
||||
const std::vector<double> dummy;
|
||||
|
||||
for (size_t i = 0; i < eps_kw.size(); ++i) {
|
||||
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
|
||||
eps_kw[i], eps_vec[i]);
|
||||
}
|
||||
|
||||
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[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,
|
||||
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,
|
||||
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,
|
||||
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,
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize saturation scaling parameter
|
||||
template<class T>
|
||||
void SaturationPropsFromDeck::initEPSKey(Opm::DeckConstPtr deck,
|
||||
Opm::EclipseStateConstPtr eclipseState,
|
||||
int number_of_cells,
|
||||
const int* global_cell,
|
||||
const T& begin_cell_centroid,
|
||||
int dimensions,
|
||||
const std::string& keyword,
|
||||
std::vector<double>& scaleparam)
|
||||
{
|
||||
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 = {
|
||||
{"SWL", 1}, {"SWCR", 2}, {"SWU", 3}, {"SGL", 4},
|
||||
{"SGCR", 5}, {"SGU", 6}, {"SOWCR", 7}, {"SOGCR", 8},
|
||||
{"ISWL", 1}, {"ISWCR", 2}, {"ISWU", 3}, {"ISGL", 4},
|
||||
{"ISGCR", 5}, {"ISGU", 6}, {"ISOWCR", 7}, {"ISOGCR", 8}};
|
||||
bool hasENPTVD = deck->hasKeyword("ENPTVD");
|
||||
bool hasENKRVD = deck->hasKeyword("ENKRVD");
|
||||
int itab = 0;
|
||||
std::vector<std::vector<double> > param_col;
|
||||
std::vector<std::vector<double> > depth_col;
|
||||
std::vector<std::string> col_names;
|
||||
std::shared_ptr<const TableManager> tables = eclipseState->getTableManager();
|
||||
|
||||
// Active keyword assigned default values for each cell (in case of possible box-wise assignment)
|
||||
if ((keyword[0] == 'S' && (useStateKeyword || hasENPTVD)) || (keyword[1] == 'S' && useStateKeyword) ) {
|
||||
if (useAqua && (useStateKeyword || columnIsMasked_(deck, "ENPTVD", kw2tab.find(keyword)->second-1))) {
|
||||
itab = kw2tab.find(keyword)->second;
|
||||
scaleparam.resize(number_of_cells);
|
||||
}
|
||||
if (!useKeyword && itab > 0) {
|
||||
const auto& enptvdTables = tables->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) {
|
||||
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) ) {
|
||||
if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) {
|
||||
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENKRVD", 0))) {
|
||||
itab = 1;
|
||||
scaleparam.resize(number_of_cells);
|
||||
for (int i=0; i<number_of_cells; ++i)
|
||||
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] = 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] = 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] = 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] = 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] = 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] = satfunc_[cell_to_func_[i]].getSatFuncBase().krorg_;
|
||||
}
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
|
||||
}
|
||||
if (!useKeyword && itab > 0) {
|
||||
const auto& enkrvdTables = tables->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) {
|
||||
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') ) {
|
||||
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] = 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] = satfunc_[cell_to_func_[i]].getSatFuncBase().pcgmax_;
|
||||
}
|
||||
}
|
||||
|
||||
if (scaleparam.empty()) {
|
||||
return;
|
||||
}
|
||||
|
||||
if (useKeyword || useStateKeyword) {
|
||||
// Keyword values from deck
|
||||
std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
|
||||
const int* gc = global_cell;
|
||||
std::vector<double> val;
|
||||
if (keyword[0] == 'S' || keyword[1] == 'S') { // Saturation from EclipseState
|
||||
val = eclipseState->getDoubleGridProperty(keyword)->getData();
|
||||
} else {
|
||||
val = deck->getKeyword(keyword)->getSIDoubleData(); //KR and PC directly from deck.
|
||||
}
|
||||
for (int c = 0; c < int(scaleparam.size()); ++c) {
|
||||
const int deck_pos = (gc == NULL) ? c : gc[c];
|
||||
scaleparam[c] = val[deck_pos];
|
||||
}
|
||||
}
|
||||
|
||||
if (itab > 0) {
|
||||
const int dim = dimensions;
|
||||
std::vector<int> endnum;
|
||||
if ( deck->hasKeyword("ENDNUM")) {
|
||||
const std::vector<int>& e =
|
||||
deck->getKeyword("ENDNUM")->getIntData();
|
||||
endnum.resize(number_of_cells);
|
||||
const int* gc = global_cell;
|
||||
for (int cell = 0; cell < number_of_cells; ++cell) {
|
||||
const int deck_pos = (gc == NULL) ? cell : gc[cell];
|
||||
endnum[cell] = e[deck_pos] - 1; // Deck value zero prevents scaling via ENPTVD/ENKRVD
|
||||
}
|
||||
}
|
||||
else {
|
||||
// Default deck value is one
|
||||
endnum.assign(number_of_cells, 0);
|
||||
}
|
||||
if (keyword[0] == 'S' || keyword[1] == 'S') { // From EclipseState
|
||||
for (int cell = 0; cell < number_of_cells; ++cell) {
|
||||
if (!std::isfinite(scaleparam[cell]) && endnum[cell] >= 0 && param_col[endnum[cell]][0] >= 0.0) {
|
||||
double zc = UgGridHelpers
|
||||
::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
|
||||
dim-1);
|
||||
if (zc >= depth_col[endnum[cell]].front() && zc <= depth_col[endnum[cell]].back()) { //don't want extrap outside depth interval
|
||||
scaleparam[cell] = linearInterpolation(depth_col[endnum[cell]], param_col[endnum[cell]], zc);
|
||||
}
|
||||
} else if (!std::isfinite(scaleparam[cell]) && endnum[cell] >= 0) {
|
||||
// As of 1/9-2014: Reflects remaining work on opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp ...
|
||||
OPM_THROW(std::runtime_error, " -- Inconsistent EclipseState: '" << keyword << "' (ENPTVD)");
|
||||
}
|
||||
}
|
||||
} else { //KR and PC from deck.
|
||||
for (int cell = 0; cell < number_of_cells; ++cell) {
|
||||
if (endnum[cell] >= 0 && param_col[endnum[cell]][0] >= 0.0) {
|
||||
double zc = UgGridHelpers
|
||||
::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
|
||||
dim-1);
|
||||
if (zc >= depth_col[endnum[cell]].front() && zc <= depth_col[endnum[cell]].back()) { //don't want extrap outside depth interval
|
||||
scaleparam[cell] = linearInterpolation(depth_col[endnum[cell]], param_col[endnum[cell]], zc);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// std::cout << keyword << ":" << std::endl;
|
||||
// for (int c = 0; c < int(scaleparam.size()); ++c) {
|
||||
// std::cout << c << " " << scaleparam[c] << std::endl;
|
||||
// }
|
||||
|
||||
}
|
||||
|
||||
// Saturation scaling
|
||||
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)
|
||||
const double scr_tab, // critical saturation
|
||||
const double su_tab, // maximum saturation (for krow/krog calculations this is minimum water/gas saturation)
|
||||
const double sxcr_tab, // second critical saturation (not used for 2pt scaling)
|
||||
const double s0_tab, // threephase complementary minimum saturation (-1.0 indicates 2-phase)
|
||||
const double krsr_tab, // relperm at displacing critical saturation
|
||||
const double krmax_tab, // relperm at maximum saturation
|
||||
const double pcmax_tab, // cap-pres at maximum saturation (zero => no scaling)
|
||||
const std::vector<double>& sl, // For krow/krog calculations this is not used
|
||||
const std::vector<double>& scr,
|
||||
const std::vector<double>& su, // For krow/krog calculations this is SWL/SGL
|
||||
const std::vector<double>& sxcr,
|
||||
const std::vector<double>& s0,
|
||||
const std::vector<double>& krsr,
|
||||
const std::vector<double>& krmax,
|
||||
const std::vector<double>& pcmax) // For krow/krog calculations this is not used
|
||||
{
|
||||
if (scr.empty() && su.empty() && (sxcr.empty() || !do_3pt_) && s0.empty()) {
|
||||
data.doNotScale = true;
|
||||
data.smin = sl_tab;
|
||||
if (oil) {
|
||||
data.smax = (s0_tab < 0.0) ? 1.0 - su_tab : 1.0 - su_tab - s0_tab;
|
||||
} else {
|
||||
data.smax = su_tab;
|
||||
}
|
||||
data.scr = scr_tab;
|
||||
} else {
|
||||
data.doNotScale = false;
|
||||
data.do_3pt = do_3pt_;
|
||||
double s_r;
|
||||
if (s0_tab < 0.0) { // 2phase
|
||||
s_r = 1.0-sxcr_tab;
|
||||
if (do_3pt_) data.sr = sxcr.empty() ? s_r : 1.0-sxcr[cell];
|
||||
} else { // 3phase
|
||||
s_r = 1.0-sxcr_tab-s0_tab;
|
||||
if (do_3pt_)data.sr = 1.0 - (sxcr.empty() ? sxcr_tab : sxcr[cell])
|
||||
- (s0.empty() ? s0_tab : s0[cell]);
|
||||
}
|
||||
data.scr = scr.empty() ? scr_tab : scr[cell];
|
||||
double s_max = su_tab;
|
||||
if (oil) {
|
||||
data.smin = sl_tab;
|
||||
if (s0_tab < 0.0) { // 2phase
|
||||
s_max = 1.0 - su_tab;
|
||||
data.smax = 1.0 - (su.empty() ? su_tab : su[cell]);
|
||||
} else { // 3phase
|
||||
s_max = 1.0 - su_tab - s0_tab;
|
||||
data.smax = 1.0 - (su.empty() ? su_tab : su[cell])
|
||||
- (s0.empty() ? s0_tab : s0[cell]);
|
||||
}
|
||||
} else {
|
||||
data.smin = sl.empty() ? sl_tab : sl[cell];
|
||||
data.smax = su.empty() ? su_tab : su[cell];
|
||||
}
|
||||
if (do_3pt_) {
|
||||
data.slope1 = (s_r-scr_tab)/(data.sr-data.scr);
|
||||
data.slope2 = (s_max-s_r)/(data.smax-data.sr);
|
||||
} else {
|
||||
data.slope2 = data.slope1 = (s_max-scr_tab)/(data.smax-data.scr);
|
||||
// Inv transform of tabulated critical displacing saturation to prepare for possible value scaling (krwr etc)
|
||||
data.sr = data.scr + (s_r-scr_tab)*(data.smax-data.scr)/(s_max-scr_tab);
|
||||
}
|
||||
}
|
||||
|
||||
data.doKrMax = !krmax.empty();
|
||||
data.doKrCrit = !krsr.empty();
|
||||
data.doSatInterp = false;
|
||||
data.krsr = krsr.empty() ? krsr_tab : krsr[cell];
|
||||
data.krmax = krmax.empty() ? krmax_tab : krmax[cell];
|
||||
data.krSlopeCrit = data.krsr/krsr_tab;
|
||||
data.krSlopeMax = data.krmax/krmax_tab;
|
||||
if (data.doKrCrit) {
|
||||
if (data.sr > data.smax-1.0e-6) {
|
||||
//Ignore krsr and do two-point (one might consider combining krsr and krmax linearly between scr and smax ... )
|
||||
data.doKrCrit = false;
|
||||
} else if (std::fabs(krmax_tab- krsr_tab) > 1.0e-6) { // interpolate kr
|
||||
data.krSlopeMax = (data.krmax-data.krsr)/(krmax_tab-krsr_tab);
|
||||
} else { // interpolate sat
|
||||
data.doSatInterp = true;
|
||||
data.krSlopeMax = (data.krmax-data.krsr)/(data.smax-data.sr);
|
||||
}
|
||||
}
|
||||
|
||||
if (std::fabs(pcmax_tab) < 1.0e-8 || pcmax.empty() || pcmax_tab*pcmax[cell] < 0.0) {
|
||||
data.pcFactor = 1.0;
|
||||
} else {
|
||||
data.pcFactor = pcmax[cell]/pcmax_tab;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED
|
||||
@@ -20,7 +20,9 @@
|
||||
#ifndef OPM_EXPLICIT_ARRAYS_FLUID_STATE_HEADER_INCLUDED
|
||||
#define OPM_EXPLICIT_ARRAYS_FLUID_STATE_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/simulator/BlackoilState.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
|
||||
#include <array>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
@@ -35,9 +37,10 @@ class ExplicitArraysFluidState
|
||||
{
|
||||
public:
|
||||
typedef double Scalar;
|
||||
enum { numPhases = BlackoilPhases::MaxNumPhases };
|
||||
|
||||
explicit ExplicitArraysFluidState(const unsigned int num_phases)
|
||||
: numPhases_(num_phases)
|
||||
explicit ExplicitArraysFluidState(const PhaseUsage& phaseUsage)
|
||||
: phaseUsage_(phaseUsage)
|
||||
{}
|
||||
|
||||
/*!
|
||||
@@ -47,7 +50,17 @@ public:
|
||||
* index.
|
||||
*/
|
||||
void setIndex(unsigned arrayIdx)
|
||||
{ arrayIdx_ = arrayIdx; }
|
||||
{
|
||||
int np = phaseUsage_.num_phases;
|
||||
for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
|
||||
if (!phaseUsage_.phase_used[phaseIdx]) {
|
||||
sats_[phaseIdx] = 0.0;
|
||||
}
|
||||
else {
|
||||
sats_[phaseIdx] = saturations_[np*arrayIdx + phaseUsage_.phase_pos[phaseIdx]];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Set the array containing the phase saturations.
|
||||
@@ -60,35 +73,18 @@ public:
|
||||
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]; }
|
||||
{ return sats_[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
|
||||
// TODO (?) temperature, pressure, composition, etc
|
||||
|
||||
private:
|
||||
const PhaseUsage phaseUsage_;
|
||||
const double* saturations_;
|
||||
const double* temperature_;
|
||||
|
||||
unsigned arrayIdx_;
|
||||
unsigned numPhases_;
|
||||
std::array<Scalar, BlackoilPhases::MaxNumPhases> sats_;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
117
opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp
Normal file
117
opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp
Normal file
@@ -0,0 +1,117 @@
|
||||
/*
|
||||
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_SAT_DERIVATIVES_FLUID_STATE_HEADER_INCLUDED
|
||||
#define OPM_EXPLICIT_ARRAYS_SAT_DERIVATIVES_FLUID_STATE_HEADER_INCLUDED
|
||||
|
||||
#include <opm/material/localad/Evaluation.hpp>
|
||||
#include <opm/material/localad/Math.hpp>
|
||||
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
|
||||
#include <array>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
// this class does not need to be implemented. its only purpose is to make compiler
|
||||
// messages for local-AD framework more explicit (because the class name of the
|
||||
// Evaluation will include a tag name so that it is clear which derivatives are handled).
|
||||
class SaturationDerivativesTag
|
||||
{};
|
||||
|
||||
/*!
|
||||
* \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::ExplicitArraysFluidState except for the fact that it
|
||||
* uses opm-material's local automatic differentiation framework to allow implicit
|
||||
* treatment of saturation derivatives.
|
||||
*/
|
||||
class ExplicitArraysSatDerivativesFluidState
|
||||
{
|
||||
public:
|
||||
enum { numPhases = BlackoilPhases::MaxNumPhases };
|
||||
enum { numComponents = 3 };
|
||||
|
||||
typedef Opm::LocalAd::Evaluation<double, SaturationDerivativesTag, numPhases> Evaluation;
|
||||
typedef Evaluation Scalar;
|
||||
|
||||
ExplicitArraysSatDerivativesFluidState(const PhaseUsage& phaseUsage)
|
||||
: phaseUsage_(phaseUsage)
|
||||
{
|
||||
globalSaturationArray_ = 0;
|
||||
|
||||
// initialize the evaluation objects for the saturations
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
|
||||
saturation_[phaseIdx] = Evaluation::createVariable(0.0, phaseIdx);
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \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)
|
||||
{
|
||||
int np = phaseUsage_.num_phases;
|
||||
|
||||
// copy the saturations values from the global value. the derivatives do not need
|
||||
// to be modified for these...
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
|
||||
if (!phaseUsage_.phase_used[phaseIdx]) {
|
||||
saturation_[phaseIdx].value = 0.0;
|
||||
}
|
||||
else {
|
||||
saturation_[phaseIdx].value = globalSaturationArray_[np*arrayIdx + phaseUsage_.phase_pos[phaseIdx]];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \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)
|
||||
{ globalSaturationArray_ = saturations; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the saturation of a phase for the current cell index.
|
||||
*/
|
||||
const Evaluation& saturation(int phaseIdx) const
|
||||
{ return saturation_[phaseIdx]; }
|
||||
|
||||
// TODO (?) temperature, pressure, composition, etc
|
||||
|
||||
private:
|
||||
const PhaseUsage phaseUsage_;
|
||||
const double* globalSaturationArray_;
|
||||
|
||||
std::array<Evaluation, numPhases> saturation_;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
@@ -693,7 +693,9 @@ namespace Opm
|
||||
// Adjust oil pressure according to gas saturation and cap pressure
|
||||
double pc[BlackoilPhases::MaxNumPhases];
|
||||
double sat[BlackoilPhases::MaxNumPhases];
|
||||
sat[waterpos] = sw;
|
||||
sat[gaspos] = sg;
|
||||
sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
|
||||
}
|
||||
@@ -704,6 +706,9 @@ namespace Opm
|
||||
double sat[BlackoilPhases::MaxNumPhases];
|
||||
double threshold_sat = 1.0e-6;
|
||||
|
||||
sat[waterpos] = smax[waterpos];
|
||||
sat[gaspos] = smax[gaspos];
|
||||
sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos];
|
||||
if (sw > smax[waterpos]-threshold_sat ) {
|
||||
sat[waterpos] = smax[waterpos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
|
||||
@@ -15,7 +15,7 @@ Source0: https://github.com/OPM/%{name}/archive/release/%{version}/%{tag}
|
||||
BuildRequires: blas-devel lapack-devel dune-common-devel
|
||||
BuildRequires: git suitesparse-devel cmake28 doxygen bc
|
||||
BuildRequires: tinyxml-devel dune-istl-devel ert.ecl-devel
|
||||
BuildRequires: opm-parser-devel boost148-devel
|
||||
BuildRequires: opm-parser-devel opm-material-devel boost148-devel
|
||||
%{?el5:BuildRequires: gcc44 gcc44-c++}
|
||||
%{!?el5:BuildRequires: gcc gcc-c++}
|
||||
BuildRoot: %{_tmppath}/%{name}-%{version}-build
|
||||
|
||||
@@ -99,10 +99,9 @@ BOOST_AUTO_TEST_CASE (GwsegStandard)
|
||||
|
||||
double krw[11] = {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7};
|
||||
double kro[11] = {1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0};
|
||||
double DkrwDsw[11] = {0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0};
|
||||
double DkroDsw[11] = {-2.0, -2.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0};
|
||||
double DkroDsg[11] = {-5.0, -5.0, -3.0, -2.0, -0.66666666666666741, -0.75, -0.8,
|
||||
-0.83333333333333237, 0.14285714285714296, 0.0, 0.0};
|
||||
double DkrwDsw[11] = {0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0};
|
||||
double DkroDsw[11] = {0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0};
|
||||
double DkroDsg[11] = {0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0};
|
||||
|
||||
const double reltol = 1.0e-6;
|
||||
for (int i=0; i<n; ++i) {
|
||||
@@ -113,7 +112,8 @@ BOOST_AUTO_TEST_CASE (GwsegStandard)
|
||||
CHECK(dkrds[i*np*np+np*gpos+opos], DkroDsg[i], reltol);
|
||||
}
|
||||
|
||||
/*
|
||||
|
||||
/*
|
||||
std::cout << std::setw(12) << "sw";
|
||||
std::cout << std::setw(12) << "so";
|
||||
std::cout << std::setw(12) << "sg";
|
||||
@@ -185,9 +185,9 @@ BOOST_AUTO_TEST_CASE (GwsegEPSBase)
|
||||
|
||||
double krw[11] = {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7};
|
||||
double kro[11] = {1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0};
|
||||
double DkrwDsw[11] = {0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0};
|
||||
double DkroDsw[11] = {-2.0, -2.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0};
|
||||
double DkroDsg[11] = {-5.0, -5.0, -3.0, -2.0,-0.66666666666666741, -0.75, -0.8, -0.83333333333333237, 0.14285714285714296, 0.0, 0.0};
|
||||
double DkrwDsw[11] = {0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0};
|
||||
double DkroDsw[11] = {0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0};
|
||||
double DkroDsg[11] = {0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0};
|
||||
|
||||
const double reltol = 1.0e-6;
|
||||
for (int i=0; i<n; ++i) {
|
||||
@@ -197,7 +197,7 @@ BOOST_AUTO_TEST_CASE (GwsegEPSBase)
|
||||
CHECK(dkrds[i*np*np+opos], DkroDsw[i], reltol);
|
||||
CHECK(dkrds[i*np*np+np*gpos+opos], DkroDsg[i], reltol);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
std::cout << std::setw(12) << "sw";
|
||||
std::cout << std::setw(12) << "so";
|
||||
@@ -278,30 +278,35 @@ BOOST_AUTO_TEST_CASE (GwsegEPS_A)
|
||||
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
|
||||
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
|
||||
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0}};
|
||||
double DkrwDsw[ncell][n] = {{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
|
||||
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
|
||||
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
|
||||
{0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0},
|
||||
{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
|
||||
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
|
||||
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
|
||||
{0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0}};
|
||||
double DkroDsw[ncell][n] = {{0.0, 0.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0},
|
||||
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
|
||||
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
|
||||
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
|
||||
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
|
||||
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
|
||||
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
|
||||
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0}};
|
||||
double DkroDsg[ncell][n] = {{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
|
||||
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
|
||||
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
|
||||
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
|
||||
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
|
||||
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
|
||||
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
|
||||
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0}};
|
||||
double DkrwDsw[ncell][n] = {
|
||||
{0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0},
|
||||
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
|
||||
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 1.4, 0},
|
||||
{0, 0, 0, 0, 0, 2.33333333333, 2.33333333333, 0, 0, 0, 0},
|
||||
{0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0},
|
||||
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
|
||||
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 1.4, 0},
|
||||
{0, 0, 0, 0, 0, 2.33333333333, 2.33333333333, 0, 0, 0, 0} };
|
||||
|
||||
double DkroDsw[ncell][n] = {
|
||||
{0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0},
|
||||
{0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0},
|
||||
{0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0},
|
||||
{0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0},
|
||||
{0, 0, 0, -2.33333333333, -2.33333333333, -1.16666666667, -1.16666666667, -1.16666666667, -1.16666666667, 0, 0},
|
||||
{0, 0, 0, -2.33333333333, -2.33333333333, -1.16666666667, -1.16666666667, -1.16666666667, -1.16666666667, 0, 0},
|
||||
{0, 0, 0, -2.33333333333, -2.33333333333, -1.16666666667, -1.16666666667, -1.16666666667, -1.16666666667, 0, 0},
|
||||
{0, 0, 0, -2.33333333333, -2.33333333333, -1.16666666667, -1.16666666667, -1.16666666667, -1.16666666667, 0, 0} };
|
||||
|
||||
double DkroDsg[ncell][n] = {
|
||||
{0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0},
|
||||
{0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0},
|
||||
{0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0},
|
||||
{0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0},
|
||||
{0, 0, 0, -3.16666666667, -2.16666666667, -0.833333333333, -0.916666666667, -0.966666666667, -1, 0, 0},
|
||||
{0, 0, 0, -3.16666666667, -2.16666666667, -0.833333333333, -0.916666666667, -0.966666666667, -1, 0, 0},
|
||||
{0, 0, 0, -3.16666666667, -2.16666666667, -0.833333333333, -0.916666666667, -0.966666666667, -1, 0, 0},
|
||||
{0, 0, 0, -3.16666666667, -2.16666666667, -0.833333333333, -0.916666666667, -0.966666666667, -1, 0, 0} };
|
||||
|
||||
for (int icell=0; icell<ncell; ++icell) {
|
||||
for (int i=0; i<n; ++i) {
|
||||
@@ -525,30 +530,36 @@ BOOST_AUTO_TEST_CASE (GwsegEPS_C)
|
||||
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
|
||||
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
|
||||
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0}};
|
||||
double DkrwDsw[ncell][n] = {{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
|
||||
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
|
||||
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
|
||||
{0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0},
|
||||
{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
|
||||
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
|
||||
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
|
||||
{0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0}};
|
||||
double DkroDsw[ncell][n] = {{0.0, 0.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0},
|
||||
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
|
||||
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
|
||||
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
|
||||
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
|
||||
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
|
||||
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
|
||||
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0}};
|
||||
double DkroDsg[ncell][n] = {{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
|
||||
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
|
||||
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
|
||||
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
|
||||
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
|
||||
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
|
||||
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
|
||||
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0}};
|
||||
double DkrwDsw[ncell][n] = {
|
||||
{0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0},
|
||||
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
|
||||
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 1.4, 0},
|
||||
{0, 0, 0, 0, 0, 2.33333333333, 2.33333333333, 0, 0, 0, 0},
|
||||
{0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0},
|
||||
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 1.4, 0, 0, 0},
|
||||
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 1.4, 0},
|
||||
{0, 0, 0, 0, 0, 2.33333333333, 2.33333333333, 0, 0, 0, 0},
|
||||
};
|
||||
double DkroDsw[ncell][n] = {
|
||||
{0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0},
|
||||
{0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0},
|
||||
{0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0},
|
||||
{0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0},
|
||||
{0, 0, 0, -2.33333333333, -2.33333333333, -1.16666666667, -1.16666666667, -1.16666666667, -1.16666666667, 0, 0},
|
||||
{0, 0, 0, -2.33333333333, -2.33333333333, -1.16666666667, -1.16666666667, -1.16666666667, -1.16666666667, 0, 0},
|
||||
{0, 0, 0, -2.33333333333, -2.33333333333, -1.16666666667, -1.16666666667, -1.16666666667, -1.16666666667, 0, 0},
|
||||
{0, 0, 0, -2.33333333333, -2.33333333333, -1.16666666667, -1.16666666667, -1.16666666667, -1.16666666667, 0, 0},
|
||||
};
|
||||
double DkroDsg[ncell][n] = {
|
||||
{0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0},
|
||||
{0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0},
|
||||
{0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0},
|
||||
{0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0},
|
||||
{0, 0, 0, -3.16666666667, -2.16666666667, -0.833333333333, -0.916666666667, -0.966666666667, -1, 0, 0},
|
||||
{0, 0, 0, -3.16666666667, -2.16666666667, -0.833333333333, -0.916666666667, -0.966666666667, -1, 0, 0},
|
||||
{0, 0, 0, -3.16666666667, -2.16666666667, -0.833333333333, -0.916666666667, -0.966666666667, -1, 0, 0},
|
||||
{0, 0, 0, -3.16666666667, -2.16666666667, -0.833333333333, -0.916666666667, -0.966666666667, -1, 0, 0},
|
||||
};
|
||||
|
||||
for (int icell=0; icell<ncell; ++icell) {
|
||||
for (int i=0; i<n; ++i) {
|
||||
@@ -568,9 +579,7 @@ BOOST_AUTO_TEST_CASE (GwsegEPS_C)
|
||||
CHECK(dkrds[i*np*np+opos], DkroDsw[icell][i], reltol);
|
||||
CHECK(dkrds[i*np*np+np*gpos+opos], DkroDsg[icell][i], reltol);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
BOOST_AUTO_TEST_CASE (GwsegEPS_D)
|
||||
@@ -617,10 +626,9 @@ BOOST_AUTO_TEST_CASE (GwsegEPS_D)
|
||||
|
||||
double krw[11] = {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7};
|
||||
double kro[11] = {1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0};
|
||||
double DkrwDsw[11] = {0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0};
|
||||
double DkroDsw[11] = {-2.0, -2.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0};
|
||||
double DkroDsg[11] = {-5.0, -5.0, -3.0, -2.0, -0.66666666666666741, -0.75, -0.8,
|
||||
-0.83333333333333237, 0.14285714285714296, 0.0, 0.0};
|
||||
double DkrwDsw[11] = {0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0};
|
||||
double DkroDsw[11] = {0, 0, -2, -2, -2, -1, -1, -1, -1, 0, 0};
|
||||
double DkroDsg[11] = {0, 0, -3, -2, -1.66666666667, -0.75, -0.8, -0.833333333333, -0.857142857143, 3.46944695195e-17, 0};
|
||||
|
||||
const double reltol = 1.0e-6;
|
||||
for (int i=0; i<n; ++i) {
|
||||
|
||||
Reference in New Issue
Block a user