Implement multi-region PVT for all property classes
since they are not using a single PVT table anymore, their "Single" prefix has been removed...
This commit is contained in:
Andreas Lauser
parent
8afa34db4d
commit
fbd8d42e8d
@ -84,11 +84,11 @@ list (APPEND MAIN_SOURCE_FILES
|
||||
opm/core/props/pvt/BlackoilPvtProperties.cpp
|
||||
opm/core/props/pvt/PvtPropertiesBasic.cpp
|
||||
opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
|
||||
opm/core/props/pvt/SinglePvtDead.cpp
|
||||
opm/core/props/pvt/SinglePvtDeadSpline.cpp
|
||||
opm/core/props/pvt/SinglePvtInterface.cpp
|
||||
opm/core/props/pvt/SinglePvtLiveGas.cpp
|
||||
opm/core/props/pvt/SinglePvtLiveOil.cpp
|
||||
opm/core/props/pvt/PvtDead.cpp
|
||||
opm/core/props/pvt/PvtDeadSpline.cpp
|
||||
opm/core/props/pvt/PvtInterface.cpp
|
||||
opm/core/props/pvt/PvtLiveGas.cpp
|
||||
opm/core/props/pvt/PvtLiveOil.cpp
|
||||
opm/core/props/rock/RockBasic.cpp
|
||||
opm/core/props/rock/RockCompressibility.cpp
|
||||
opm/core/props/rock/RockFromDeck.cpp
|
||||
@ -311,12 +311,12 @@ list (APPEND PUBLIC_HEADER_FILES
|
||||
opm/core/props/pvt/BlackoilPvtProperties.hpp
|
||||
opm/core/props/pvt/PvtPropertiesBasic.hpp
|
||||
opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp
|
||||
opm/core/props/pvt/SinglePvtConstCompr.hpp
|
||||
opm/core/props/pvt/SinglePvtDead.hpp
|
||||
opm/core/props/pvt/SinglePvtDeadSpline.hpp
|
||||
opm/core/props/pvt/SinglePvtInterface.hpp
|
||||
opm/core/props/pvt/SinglePvtLiveGas.hpp
|
||||
opm/core/props/pvt/SinglePvtLiveOil.hpp
|
||||
opm/core/props/pvt/PvtConstCompr.hpp
|
||||
opm/core/props/pvt/PvtDead.hpp
|
||||
opm/core/props/pvt/PvtDeadSpline.hpp
|
||||
opm/core/props/pvt/PvtInterface.hpp
|
||||
opm/core/props/pvt/PvtLiveGas.hpp
|
||||
opm/core/props/pvt/PvtLiveOil.hpp
|
||||
opm/core/props/rock/RockBasic.hpp
|
||||
opm/core/props/rock/RockCompressibility.hpp
|
||||
opm/core/props/rock/RockFromDeck.hpp
|
||||
|
||||
@ -23,8 +23,14 @@ namespace Opm
|
||||
const parameter::ParameterGroup& param,
|
||||
bool init_rock)
|
||||
{
|
||||
init(deck, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
|
||||
param, init_rock);
|
||||
init(deck,
|
||||
number_of_cells,
|
||||
global_cell,
|
||||
cart_dims,
|
||||
begin_cell_centroids,
|
||||
dimension,
|
||||
param,
|
||||
init_rock);
|
||||
}
|
||||
|
||||
template<class CentroidIterator>
|
||||
@ -39,7 +45,7 @@ namespace Opm
|
||||
if (init_rock){
|
||||
rock_.init(deck, number_of_cells, global_cell, cart_dims);
|
||||
}
|
||||
pvt_.init(deck, /*numSamples=*/0);
|
||||
pvt_.init(deck, /*numSamples=*/0, number_of_cells, global_cell);
|
||||
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
|
||||
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
|
||||
satprops_.reset(ptr);
|
||||
@ -68,7 +74,7 @@ namespace Opm
|
||||
}
|
||||
|
||||
const int pvt_samples = param.getDefault("pvt_tab_size", 200);
|
||||
pvt_.init(deck, pvt_samples);
|
||||
pvt_.init(deck, pvt_samples, number_of_cells, global_cell);
|
||||
|
||||
// Unfortunate lack of pointer smartness here...
|
||||
const int sat_samples = param.getDefault("sat_tab_size", 200);
|
||||
|
||||
@ -21,11 +21,11 @@
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtDead.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtConstCompr.hpp>
|
||||
#include <opm/core/props/pvt/PvtDead.hpp>
|
||||
#include <opm/core/props/pvt/PvtDeadSpline.hpp>
|
||||
#include <opm/core/props/pvt/PvtLiveOil.hpp>
|
||||
#include <opm/core/props/pvt/PvtLiveGas.hpp>
|
||||
#include <opm/core/props/pvt/PvtConstCompr.hpp>
|
||||
#include <opm/core/props/phaseUsageFromDeck.hpp>
|
||||
#include <opm/core/utility/Units.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
@ -42,93 +42,107 @@ namespace Opm
|
||||
{
|
||||
}
|
||||
|
||||
void BlackoilPvtProperties::init(Opm::DeckConstPtr deck, int samples)
|
||||
void BlackoilPvtProperties::init(Opm::DeckConstPtr deck,
|
||||
int numSamples,
|
||||
int numCompressedCells,
|
||||
const int *compressedToCartesianCellIdx)
|
||||
{
|
||||
// If we need multiple regions, this class and the SinglePvt* classes must change.
|
||||
region_number_ = 0;
|
||||
|
||||
phase_usage_ = phaseUsageFromDeck(deck);
|
||||
|
||||
// Surface densities. Accounting for different orders in eclipse and our code.
|
||||
if (deck->hasKeyword("DENSITY")) {
|
||||
Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY");
|
||||
Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY");
|
||||
int numRegions = densityKeyword->size();
|
||||
|
||||
densities_.resize(numRegions);
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
if (phase_usage_.phase_used[Liquid]) {
|
||||
densities_[phase_usage_.phase_pos[Liquid]]
|
||||
= densityKeyword->getRecord(region_number_)->getItem("OIL")->getSIDouble(0);
|
||||
densities_[regionIdx][phase_usage_.phase_pos[Liquid]]
|
||||
= densityKeyword->getRecord(regionIdx)->getItem("OIL")->getSIDouble(0);
|
||||
}
|
||||
if (phase_usage_.phase_used[Aqua]) {
|
||||
densities_[phase_usage_.phase_pos[Aqua]]
|
||||
= densityKeyword->getRecord(region_number_)->getItem("WATER")->getSIDouble(0);
|
||||
densities_[regionIdx][phase_usage_.phase_pos[Aqua]]
|
||||
= densityKeyword->getRecord(regionIdx)->getItem("WATER")->getSIDouble(0);
|
||||
}
|
||||
if (phase_usage_.phase_used[Vapour]) {
|
||||
densities_[phase_usage_.phase_pos[Vapour]]
|
||||
= densityKeyword->getRecord(region_number_)->getItem("GAS")->getSIDouble(0);
|
||||
densities_[regionIdx][phase_usage_.phase_pos[Vapour]]
|
||||
= densityKeyword->getRecord(regionIdx)->getItem("GAS")->getSIDouble(0);
|
||||
}
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");
|
||||
}
|
||||
|
||||
// Set the properties.
|
||||
// first, calculate the PVT table index for each compressed
|
||||
// cell. This array is required to construct the PVT classes
|
||||
// below.
|
||||
Opm::extractPvtTableIndex(pvtTableIdx_,
|
||||
deck,
|
||||
numCompressedCells,
|
||||
compressedToCartesianCellIdx);
|
||||
|
||||
// Resize the property objects container
|
||||
props_.resize(phase_usage_.num_phases);
|
||||
|
||||
// Water PVT
|
||||
if (phase_usage_.phase_used[Aqua]) {
|
||||
if (deck->hasKeyword("PVTW")) {
|
||||
Opm::PvtwTable pvtwTable(deck->getKeyword("PVTW"));
|
||||
// if water is used, we require the presence of the "PVTW"
|
||||
// keyword for now...
|
||||
std::shared_ptr<PvtConstCompr> pvtw(new PvtConstCompr);
|
||||
pvtw->initFromWater(deck->getKeyword("PVTW"), pvtTableIdx_);
|
||||
|
||||
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable));
|
||||
} else {
|
||||
// Eclipse 100 default.
|
||||
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
|
||||
}
|
||||
props_[phase_usage_.phase_pos[Aqua]] = pvtw;
|
||||
}
|
||||
// Oil PVT
|
||||
if (phase_usage_.phase_used[Liquid]) {
|
||||
// for oil, we support the "PVDO", "PVTO" and "PVCDO"
|
||||
// keywords...
|
||||
if (deck->hasKeyword("PVDO")) {
|
||||
Opm::PvdoTable pvdoTable(deck->getKeyword("PVDO"), region_number_);
|
||||
if (samples > 0) {
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
|
||||
Opm::DeckKeywordConstPtr pvdoKeyword = deck->getKeyword("PVDO");
|
||||
if (numSamples > 0) {
|
||||
auto splinePvt = std::shared_ptr<PvtDeadSpline>(new PvtDeadSpline);
|
||||
splinePvt->initFromOil(pvdoKeyword, pvtTableIdx_, numSamples);
|
||||
props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
|
||||
} else {
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable));
|
||||
auto deadPvt = std::shared_ptr<PvtDead>(new PvtDead);
|
||||
deadPvt->initFromOil(pvdoKeyword, pvtTableIdx_);
|
||||
props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
|
||||
}
|
||||
} else if (deck->hasKeyword("PVTO")) {
|
||||
Opm::PvtoTable pvtoTable(deck->getKeyword("PVTO"), /*tableIdx=*/0);
|
||||
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable));
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"), pvtTableIdx_));
|
||||
} else if (deck->hasKeyword("PVCDO")) {
|
||||
Opm::PvcdoTable pvcdoTable(deck->getKeyword("PVCDO"));
|
||||
std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
|
||||
pvcdo->initFromOil(deck->getKeyword("PVCDO"), pvtTableIdx_);
|
||||
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable));
|
||||
props_[phase_usage_.phase_pos[Liquid]] = pvcdo;
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
|
||||
OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n");
|
||||
}
|
||||
}
|
||||
// Gas PVT
|
||||
if (phase_usage_.phase_used[Vapour]) {
|
||||
// gas can be specified using the "PVDG" or "PVTG" keywords...
|
||||
if (deck->hasKeyword("PVDG")) {
|
||||
Opm::PvdgTable pvdgTable(deck->getKeyword("PVDG"), region_number_);
|
||||
if (samples > 0) {
|
||||
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
|
||||
Opm::DeckKeywordConstPtr pvdgKeyword = deck->getKeyword("PVDG");
|
||||
|
||||
if (numSamples > 0) {
|
||||
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
|
||||
splinePvt->initFromGas(pvdgKeyword, pvtTableIdx_, numSamples);
|
||||
|
||||
props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
|
||||
} else {
|
||||
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable));
|
||||
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
|
||||
deadPvt->initFromGas(pvdgKeyword, pvtTableIdx_);
|
||||
|
||||
props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
|
||||
}
|
||||
} else if (deck->hasKeyword("PVTG")) {
|
||||
Opm::PvtgTable pvtgTable(deck->getKeyword("PVTG"), /*tableIdx=*/0);
|
||||
|
||||
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable));
|
||||
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG"), pvtTableIdx_));
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
|
||||
}
|
||||
}
|
||||
|
||||
// Must inform pvt property objects of phase structure.
|
||||
for (int i = 0; i < phase_usage_.num_phases; ++i) {
|
||||
props_[i]->setPhaseConfiguration(phase_usage_.num_phases, phase_usage_.phase_pos);
|
||||
}
|
||||
}
|
||||
|
||||
const double* BlackoilPvtProperties::surfaceDensities() const
|
||||
const double* BlackoilPvtProperties::surfaceDensities(int regionIdx) const
|
||||
{
|
||||
return densities_;
|
||||
return &densities_[regionIdx][0];
|
||||
}
|
||||
|
||||
|
||||
@ -234,5 +248,4 @@ namespace Opm
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
@ -20,10 +20,11 @@
|
||||
#ifndef OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
|
||||
#define OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
|
||||
#include <opm/core/props/pvt/PvtInterface.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
|
||||
|
||||
#include <string>
|
||||
#include <memory>
|
||||
@ -37,9 +38,9 @@ namespace Opm
|
||||
/// Output arrays shall be of size n*num_phases, and must be valid
|
||||
/// before calling the method.
|
||||
/// NOTE: The difference between this interface and the one defined
|
||||
/// by SinglePvtInterface is that this collects all phases' properties,
|
||||
/// by PvtInterface is that this collects all phases' properties,
|
||||
/// and therefore the output arrays are of size n*num_phases as opposed
|
||||
/// to size n in SinglePvtInterface.
|
||||
/// to size n in PvtInterface.
|
||||
class BlackoilPvtProperties : public BlackoilPhases
|
||||
{
|
||||
public:
|
||||
@ -47,8 +48,17 @@ namespace Opm
|
||||
BlackoilPvtProperties();
|
||||
|
||||
/// Initialize from deck.
|
||||
///
|
||||
/// This method needs the mapping for compressed to cartesian
|
||||
/// cell indices because the methods which do the actual work
|
||||
/// are specified for compressed cells, but the eclipse input
|
||||
/// data is specified on cartesian cell indices...
|
||||
///
|
||||
/// \param deck An input deck from the opm-parser module.
|
||||
void init(Opm::DeckConstPtr deck, int samples);
|
||||
void init(Opm::DeckConstPtr deck,
|
||||
int samples,
|
||||
int numCompressedCells,
|
||||
const int *compressedToCartesianCellIdx);
|
||||
|
||||
/// \return Object describing the active phases.
|
||||
PhaseUsage phaseUsage() const;
|
||||
@ -68,7 +78,7 @@ namespace Opm
|
||||
|
||||
/// Densities of stock components at surface conditions.
|
||||
/// \return Array of size numPhases().
|
||||
const double* surfaceDensities() const;
|
||||
const double* surfaceDensities(int regionIdx = 0) const;
|
||||
|
||||
/// Viscosity as a function of p and z.
|
||||
void mu(const int n,
|
||||
@ -109,11 +119,19 @@ namespace Opm
|
||||
|
||||
PhaseUsage phase_usage_;
|
||||
|
||||
int region_number_;
|
||||
// The PVT properties. One object per active fluid phase.
|
||||
std::vector<std::shared_ptr<PvtInterface> > props_;
|
||||
|
||||
std::vector<std::shared_ptr<SinglePvtInterface> > props_;
|
||||
// The index of the PVT table which ought to be used for each
|
||||
// cell. Eclipse does not seem to allow specifying fluid-phase
|
||||
// specific table indices, so for the sake of simplicity, we
|
||||
// don't do this either. (if it turns out that Eclipes does in
|
||||
// fact support it or if it by some miracle gains this feature
|
||||
// in the future, this attribute needs to become a vector of
|
||||
// vectors of ints.)
|
||||
std::vector<int> pvtTableIdx_;
|
||||
|
||||
double densities_[MaxNumPhases];
|
||||
std::vector<std::array<double, MaxNumPhases> > densities_;
|
||||
mutable std::vector<double> data1_;
|
||||
mutable std::vector<double> data2_;
|
||||
};
|
||||
|
||||
291
opm/core/props/pvt/PvtConstCompr.hpp
Normal file
291
opm/core/props/pvt/PvtConstCompr.hpp
Normal file
@ -0,0 +1,291 @@
|
||||
/*
|
||||
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_PVTCONSTCOMPR_HEADER_INCLUDED
|
||||
#define OPM_PVTCONSTCOMPR_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/props/pvt/PvtInterface.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
|
||||
#include <opm/parser/eclipse/Utility/PvcdoTable.hpp>
|
||||
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
/// Class for constant compressible phases (PVTW or PVCDO).
|
||||
/// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
|
||||
/// or pressure (p) and gas resolution factor (r).
|
||||
/// For all the virtual methods, the following apply: p, r and z
|
||||
/// are expected to be of size n, size n and n*num_phases, respectively.
|
||||
/// Output arrays shall be of size n, and must be valid before
|
||||
/// calling the method.
|
||||
class PvtConstCompr : public PvtInterface
|
||||
{
|
||||
public:
|
||||
PvtConstCompr()
|
||||
{}
|
||||
|
||||
void initFromWater(Opm::DeckKeywordConstPtr pvtwKeyword, const std::vector<int> &pvtTableIdx)
|
||||
{
|
||||
pvtTableIdx_ = pvtTableIdx;
|
||||
|
||||
int numRegions = pvtwKeyword->size();
|
||||
|
||||
ref_press_.resize(numRegions);
|
||||
ref_B_.resize(numRegions);
|
||||
comp_.resize(numRegions);
|
||||
viscosity_.resize(numRegions);
|
||||
visc_comp_.resize(numRegions);
|
||||
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
|
||||
Opm::DeckRecordConstPtr pvtwRecord = pvtwKeyword->getRecord(regionIdx);
|
||||
|
||||
ref_press_[regionIdx] = pvtwRecord->getItem("P_REF")->getSIDouble(0);
|
||||
ref_B_[regionIdx] = pvtwRecord->getItem("WATER_VOL_FACTOR")->getSIDouble(0);
|
||||
comp_[regionIdx] = pvtwRecord->getItem("WATER_COMPRESSIBILITY")->getSIDouble(0);
|
||||
viscosity_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSITY")->getSIDouble(0);
|
||||
visc_comp_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSIBILITY")->getSIDouble(0);
|
||||
}
|
||||
}
|
||||
|
||||
void initFromOil(Opm::DeckKeywordConstPtr pvcdoKeyword, const std::vector<int> &pvtTableIdx)
|
||||
{
|
||||
pvtTableIdx_ = pvtTableIdx;
|
||||
|
||||
int numRegions = pvcdoKeyword->size();
|
||||
|
||||
ref_press_.resize(numRegions);
|
||||
ref_B_.resize(numRegions);
|
||||
comp_.resize(numRegions);
|
||||
viscosity_.resize(numRegions);
|
||||
visc_comp_.resize(numRegions);
|
||||
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
|
||||
Opm::DeckRecordConstPtr pvcdoRecord = pvcdoKeyword->getRecord(regionIdx);
|
||||
|
||||
ref_press_[regionIdx] = pvcdoRecord->getItem("P_REF")->getSIDouble(0);
|
||||
ref_B_[regionIdx] = pvcdoRecord->getItem("BO_REF")->getSIDouble(0);
|
||||
comp_[regionIdx] = pvcdoRecord->getItem("C_OIL_REF")->getSIDouble(0);
|
||||
viscosity_[regionIdx] = pvcdoRecord->getItem("MUO_REF")->getSIDouble(0);
|
||||
visc_comp_[regionIdx] = pvcdoRecord->getItem("OIL_VISCOSIBILITY")->getSIDouble(0);
|
||||
}
|
||||
}
|
||||
|
||||
PvtConstCompr(double visc)
|
||||
: ref_press_(1, 0.0),
|
||||
ref_B_(1, 1.0),
|
||||
comp_(1, 0.0),
|
||||
viscosity_(1, visc),
|
||||
visc_comp_(1, 0.0)
|
||||
{
|
||||
}
|
||||
|
||||
virtual ~PvtConstCompr()
|
||||
{
|
||||
}
|
||||
|
||||
virtual void mu(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_mu) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
int tableIdx = getTableIndex_(i);
|
||||
double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
|
||||
output_mu[i] = viscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
|
||||
}
|
||||
}
|
||||
|
||||
virtual void mu(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
int tableIdx = getTableIndex_(i);
|
||||
double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
|
||||
double d = (1.0 + x + 0.5*x*x);
|
||||
output_mu[i] = viscosity_[tableIdx]/d;
|
||||
output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx];
|
||||
}
|
||||
std::fill(output_dmudr, output_dmudr + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void mu(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
const PhasePresence* /*cond*/,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
int tableIdx = getTableIndex_(i);
|
||||
double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
|
||||
double d = (1.0 + x + 0.5*x*x);
|
||||
output_mu[i] = viscosity_[tableIdx]/d;
|
||||
output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx];
|
||||
}
|
||||
std::fill(output_dmudr, output_dmudr + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void B(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_B) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
int tableIdx = getTableIndex_(i);
|
||||
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
|
||||
output_B[i] = ref_B_[tableIdx]/(1.0 + x + 0.5*x*x);
|
||||
}
|
||||
}
|
||||
|
||||
virtual void dBdp(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_B,
|
||||
double* output_dBdp) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int tableIdx = getTableIndex_(i);
|
||||
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
|
||||
double d = (1.0 + x + 0.5*x*x);
|
||||
output_B[i] = ref_B_[tableIdx]/d;
|
||||
output_dBdp[i] = (-ref_B_[tableIdx]/(d*d))*(1 + x) * comp_[tableIdx];
|
||||
}
|
||||
}
|
||||
|
||||
virtual void b(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
int tableIdx = getTableIndex_(i);
|
||||
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
|
||||
double d = (1.0 + x + 0.5*x*x);
|
||||
|
||||
// b = 1/B = d/ref_B_B;
|
||||
output_b[i] = d/ref_B_[tableIdx];
|
||||
output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx];
|
||||
}
|
||||
|
||||
std::fill(output_dbdr, output_dbdr + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void b(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
const PhasePresence* /*cond*/,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
int tableIdx = getTableIndex_(i);
|
||||
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
|
||||
double d = (1.0 + x + 0.5*x*x);
|
||||
|
||||
// b = 1/B = d/ref_B_[tableIdx]B;
|
||||
output_b[i] = d/ref_B_[tableIdx];
|
||||
output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx];
|
||||
}
|
||||
std::fill(output_dbdr, output_dbdr + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void rsSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rsSat,
|
||||
double* output_drsSatdp) const
|
||||
{
|
||||
std::fill(output_rsSat, output_rsSat + n, 0.0);
|
||||
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void rvSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rvSat,
|
||||
double* output_drvSatdp) const
|
||||
{
|
||||
std::fill(output_rvSat, output_rvSat + n, 0.0);
|
||||
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void R(const int n,
|
||||
const double* /*p*/,
|
||||
const double* /*z*/,
|
||||
double* output_R) const
|
||||
{
|
||||
std::fill(output_R, output_R + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void dRdp(const int n,
|
||||
const double* /*p*/,
|
||||
const double* /*z*/,
|
||||
double* output_R,
|
||||
double* output_dRdp) const
|
||||
{
|
||||
std::fill(output_R, output_R + n, 0.0);
|
||||
std::fill(output_dRdp, output_dRdp + n, 0.0);
|
||||
}
|
||||
|
||||
private:
|
||||
int getTableIndex_(int cellIdx) const
|
||||
{
|
||||
if (pvtTableIdx_.empty())
|
||||
return 0;
|
||||
return pvtTableIdx_[cellIdx];
|
||||
}
|
||||
|
||||
std::vector<int> pvtTableIdx_;
|
||||
std::vector<double> ref_press_;
|
||||
std::vector<double> ref_B_;
|
||||
std::vector<double> comp_;
|
||||
std::vector<double> viscosity_;
|
||||
std::vector<double> visc_comp_;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif // OPM_PVTCONSTCOMPR_HEADER_INCLUDED
|
||||
|
||||
@ -19,7 +19,7 @@
|
||||
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/props/pvt/SinglePvtDead.hpp>
|
||||
#include <opm/core/props/pvt/PvtDead.hpp>
|
||||
#include <algorithm>
|
||||
|
||||
// Extra includes for debug dumping of tables.
|
||||
@ -34,72 +34,84 @@ namespace Opm
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
/// Constructor
|
||||
SinglePvtDead::SinglePvtDead(const Opm::PvdoTable& pvdoTable)
|
||||
void PvtDead::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
|
||||
const std::vector<int> &pvtTableIdx)
|
||||
{
|
||||
// Copy data
|
||||
const std::vector<double>& press = pvdoTable.getPressureColumn();
|
||||
const std::vector<double>& b = pvdoTable.getFormationFactorColumn();
|
||||
const std::vector<double>& visc = pvdoTable.getViscosityColumn();
|
||||
pvtTableIdx_ = pvtTableIdx;
|
||||
|
||||
const int sz = b.size();
|
||||
std::vector<double> bInv(sz);
|
||||
for (int i = 0; i < sz; ++i) {
|
||||
bInv[i] = 1.0 / b[i];
|
||||
int numRegions = Opm::PvdoTable::numTables(pvdoKeyword);
|
||||
|
||||
// resize the attributes of the object
|
||||
b_.resize(numRegions);
|
||||
viscosity_.resize(numRegions);
|
||||
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx);
|
||||
// Copy data
|
||||
const std::vector<double>& press = pvdoTable.getPressureColumn();
|
||||
const std::vector<double>& b = pvdoTable.getFormationFactorColumn();
|
||||
const std::vector<double>& visc = pvdoTable.getViscosityColumn();
|
||||
|
||||
const int sz = b.size();
|
||||
std::vector<double> bInv(sz);
|
||||
for (int i = 0; i < sz; ++i) {
|
||||
bInv[i] = 1.0 / b[i];
|
||||
}
|
||||
b_[regionIdx] = NonuniformTableLinear<double>(press, bInv);
|
||||
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
|
||||
}
|
||||
b_ = NonuniformTableLinear<double>(press, bInv);
|
||||
viscosity_ = NonuniformTableLinear<double>(press, visc);
|
||||
|
||||
// Dumping the created tables.
|
||||
// static int count = 0;
|
||||
// std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
|
||||
// os.precision(15);
|
||||
// os << "1/B\n\n" << one_over_B_
|
||||
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
|
||||
}
|
||||
|
||||
/// Constructor
|
||||
SinglePvtDead::SinglePvtDead(const Opm::PvdgTable& pvdgTable)
|
||||
|
||||
void PvtDead::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
|
||||
const std::vector<int> &pvtTableIdx)
|
||||
{
|
||||
// Copy data
|
||||
const std::vector<double>& press = pvdgTable.getPressureColumn();
|
||||
const std::vector<double>& b = pvdgTable.getFormationFactorColumn();
|
||||
const std::vector<double>& visc = pvdgTable.getViscosityColumn();
|
||||
pvtTableIdx_ = pvtTableIdx;
|
||||
|
||||
const int sz = b.size();
|
||||
std::vector<double> bInv(sz);
|
||||
for (int i = 0; i < sz; ++i) {
|
||||
bInv[i] = 1.0 / b[i];
|
||||
int numRegions = Opm::PvdgTable::numTables(pvdgKeyword);
|
||||
|
||||
// resize the attributes of the object
|
||||
b_.resize(numRegions);
|
||||
viscosity_.resize(numRegions);
|
||||
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx);
|
||||
|
||||
// Copy data
|
||||
const std::vector<double>& press = pvdgTable.getPressureColumn();
|
||||
const std::vector<double>& b = pvdgTable.getFormationFactorColumn();
|
||||
const std::vector<double>& visc = pvdgTable.getViscosityColumn();
|
||||
|
||||
const int sz = b.size();
|
||||
std::vector<double> bInv(sz);
|
||||
for (int i = 0; i < sz; ++i) {
|
||||
bInv[i] = 1.0 / b[i];
|
||||
}
|
||||
b_[regionIdx] = NonuniformTableLinear<double>(press, bInv);
|
||||
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
|
||||
}
|
||||
b_ = NonuniformTableLinear<double>(press, bInv);
|
||||
viscosity_ = NonuniformTableLinear<double>(press, visc);
|
||||
|
||||
// Dumping the created tables.
|
||||
// static int count = 0;
|
||||
// std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
|
||||
// os.precision(15);
|
||||
// os << "1/B\n\n" << one_over_B_
|
||||
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
|
||||
}
|
||||
|
||||
// Destructor
|
||||
SinglePvtDead::~SinglePvtDead()
|
||||
PvtDead::~PvtDead()
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
|
||||
void SinglePvtDead::mu(const int n,
|
||||
void PvtDead::mu(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_mu) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_mu[i] = viscosity_(p[i]);
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
output_mu[i] = viscosity_[regionIdx](p[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtDead::mu(const int n,
|
||||
void PvtDead::mu(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
double* output_mu,
|
||||
@ -108,14 +120,15 @@ namespace Opm
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_mu[i] = viscosity_(p[i]);
|
||||
output_dmudp[i] = viscosity_.derivative(p[i]);
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
output_mu[i] = viscosity_[regionIdx](p[i]);
|
||||
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
|
||||
}
|
||||
std::fill(output_dmudr, output_dmudr + n, 0.0);
|
||||
|
||||
}
|
||||
|
||||
void SinglePvtDead::mu(const int n,
|
||||
void PvtDead::mu(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
const PhasePresence* /*cond*/,
|
||||
@ -125,14 +138,15 @@ namespace Opm
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_mu[i] = viscosity_(p[i]);
|
||||
output_dmudp[i] = viscosity_.derivative(p[i]);
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
output_mu[i] = viscosity_[regionIdx](p[i]);
|
||||
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
|
||||
}
|
||||
std::fill(output_dmudr, output_dmudr + n, 0.0);
|
||||
|
||||
}
|
||||
|
||||
void SinglePvtDead::B(const int n,
|
||||
void PvtDead::B(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_B) const
|
||||
@ -140,11 +154,12 @@ namespace Opm
|
||||
// #pragma omp parallel for
|
||||
// B = 1/b
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_B[i] = 1.0/b_(p[i]);
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
output_B[i] = 1.0/b_[regionIdx](p[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtDead::dBdp(const int n,
|
||||
void PvtDead::dBdp(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_B,
|
||||
@ -153,12 +168,13 @@ namespace Opm
|
||||
B(n, p, 0, output_B);
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
double Bg = output_B[i];
|
||||
output_dBdp[i] = -Bg*Bg*b_.derivative(p[i]);
|
||||
output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtDead::b(const int n,
|
||||
void PvtDead::b(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
double* output_b,
|
||||
@ -168,15 +184,17 @@ namespace Opm
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_b[i] = b_(p[i]);
|
||||
output_dbdp[i] = b_.derivative(p[i]);
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
|
||||
output_b[i] = b_[regionIdx](p[i]);
|
||||
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
|
||||
|
||||
}
|
||||
std::fill(output_dbdr, output_dbdr + n, 0.0);
|
||||
|
||||
}
|
||||
|
||||
void SinglePvtDead::b(const int n,
|
||||
void PvtDead::b(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
const PhasePresence* /*cond*/,
|
||||
@ -187,15 +205,17 @@ namespace Opm
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_b[i] = b_(p[i]);
|
||||
output_dbdp[i] = b_.derivative(p[i]);
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
|
||||
output_b[i] = b_[regionIdx](p[i]);
|
||||
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
|
||||
|
||||
}
|
||||
std::fill(output_dbdr, output_dbdr + n, 0.0);
|
||||
|
||||
}
|
||||
|
||||
void SinglePvtDead::rsSat(const int n,
|
||||
void PvtDead::rsSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rsSat,
|
||||
double* output_drsSatdp) const
|
||||
@ -204,7 +224,7 @@ namespace Opm
|
||||
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
void SinglePvtDead::rvSat(const int n,
|
||||
void PvtDead::rvSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rvSat,
|
||||
double* output_drvSatdp) const
|
||||
@ -213,7 +233,7 @@ namespace Opm
|
||||
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
void SinglePvtDead::R(const int n,
|
||||
void PvtDead::R(const int n,
|
||||
const double* /*p*/,
|
||||
const double* /*z*/,
|
||||
double* output_R) const
|
||||
@ -221,7 +241,7 @@ namespace Opm
|
||||
std::fill(output_R, output_R + n, 0.0);
|
||||
}
|
||||
|
||||
void SinglePvtDead::dRdp(const int n,
|
||||
void PvtDead::dRdp(const int n,
|
||||
const double* /*p*/,
|
||||
const double* /*z*/,
|
||||
double* output_R,
|
||||
@ -17,11 +17,11 @@
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_SINGLEPVTDEAD_HEADER_INCLUDED
|
||||
#define OPM_SINGLEPVTDEAD_HEADER_INCLUDED
|
||||
#ifndef OPM_PVTDEAD_HEADER_INCLUDED
|
||||
#define OPM_PVTDEAD_HEADER_INCLUDED
|
||||
|
||||
|
||||
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
|
||||
#include <opm/core/props/pvt/PvtInterface.hpp>
|
||||
#include <opm/core/utility/NonuniformTableLinear.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
|
||||
@ -39,12 +39,16 @@ namespace Opm
|
||||
/// are expected to be of size n, size n and n*num_phases, respectively.
|
||||
/// Output arrays shall be of size n, and must be valid before
|
||||
/// calling the method.
|
||||
class SinglePvtDead : public SinglePvtInterface
|
||||
class PvtDead : public PvtInterface
|
||||
{
|
||||
public:
|
||||
SinglePvtDead(const Opm::PvdoTable& pvdoTable);
|
||||
SinglePvtDead(const Opm::PvdgTable& pvdgTable);
|
||||
virtual ~SinglePvtDead();
|
||||
PvtDead() {};
|
||||
|
||||
void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
|
||||
const std::vector<int> &pvtTableIdx);
|
||||
void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
|
||||
const std::vector<int> &pvtTableIdx);
|
||||
virtual ~PvtDead();
|
||||
|
||||
/// Viscosity as a function of p and z.
|
||||
virtual void mu(const int n,
|
||||
@ -130,11 +134,11 @@ namespace Opm
|
||||
double* output_dRdp) const;
|
||||
private:
|
||||
// PVT properties of dry gas or dead oil
|
||||
NonuniformTableLinear<double> b_;
|
||||
NonuniformTableLinear<double> viscosity_;
|
||||
std::vector<int> pvtTableIdx_;
|
||||
std::vector<NonuniformTableLinear<double> > b_;
|
||||
std::vector<NonuniformTableLinear<double> > viscosity_;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif // OPM_SINGLEPVTDEAD_HEADER_INCLUDED
|
||||
#endif // OPM_PVTDEAD_HEADER_INCLUDED
|
||||
255
opm/core/props/pvt/PvtDeadSpline.cpp
Normal file
255
opm/core/props/pvt/PvtDeadSpline.cpp
Normal file
@ -0,0 +1,255 @@
|
||||
/*
|
||||
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/props/pvt/PvtDeadSpline.hpp>
|
||||
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/SingleRecordTable.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
// Extra includes for debug dumping of tables.
|
||||
// #include <boost/lexical_cast.hpp>
|
||||
// #include <string>
|
||||
// #include <fstream>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
//------------------------------------------------------------------------
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
|
||||
PvtDeadSpline::PvtDeadSpline()
|
||||
{}
|
||||
|
||||
void PvtDeadSpline::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
|
||||
const std::vector<int> &pvtTableIdx,
|
||||
int numSamples)
|
||||
{
|
||||
pvtTableIdx_ = pvtTableIdx;
|
||||
|
||||
int numRegions = Opm::PvdoTable::numTables(pvdoKeyword);
|
||||
|
||||
// resize the attributes of the object
|
||||
b_.resize(numRegions);
|
||||
viscosity_.resize(numRegions);
|
||||
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx);
|
||||
|
||||
int numRows = pvdoTable.numRows();
|
||||
|
||||
// Copy data
|
||||
const std::vector<double> &press = pvdoTable.getPressureColumn();
|
||||
const std::vector<double> &B = pvdoTable.getFormationFactorColumn();
|
||||
const std::vector<double> &visc = pvdoTable.getViscosityColumn();
|
||||
|
||||
std::vector<double> B_inv(numRows);
|
||||
for (int i = 0; i < numRows; ++i) {
|
||||
B_inv[i] = 1.0 / B[i];
|
||||
}
|
||||
|
||||
buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]);
|
||||
buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]);
|
||||
}
|
||||
}
|
||||
|
||||
void PvtDeadSpline::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
|
||||
const std::vector<int> &pvtTableIdx,
|
||||
int numSamples)
|
||||
{
|
||||
pvtTableIdx_ = pvtTableIdx;
|
||||
|
||||
int numRegions = Opm::PvdgTable::numTables(pvdgKeyword);
|
||||
|
||||
// resize the attributes of the object
|
||||
b_.resize(numRegions);
|
||||
viscosity_.resize(numRegions);
|
||||
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx);
|
||||
|
||||
int numRows = pvdgTable.numRows();
|
||||
|
||||
// Copy data
|
||||
const std::vector<double> &press = pvdgTable.getPressureColumn();
|
||||
const std::vector<double> &B = pvdgTable.getFormationFactorColumn();
|
||||
const std::vector<double> &visc = pvdgTable.getViscosityColumn();
|
||||
|
||||
std::vector<double> B_inv(numRows);
|
||||
for (int i = 0; i < numRows; ++i) {
|
||||
B_inv[i] = 1.0 / B[i];
|
||||
}
|
||||
|
||||
buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]);
|
||||
buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]);
|
||||
}
|
||||
}
|
||||
|
||||
// Destructor
|
||||
PvtDeadSpline::~PvtDeadSpline()
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
|
||||
void PvtDeadSpline::mu(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_mu) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
output_mu[i] = viscosity_[regionIdx](p[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void PvtDeadSpline::mu(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
output_mu[i] = viscosity_[regionIdx](p[i]);
|
||||
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
|
||||
}
|
||||
std::fill(output_dmudr, output_dmudr + n, 0.0);
|
||||
}
|
||||
|
||||
void PvtDeadSpline::mu(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
const PhasePresence* /*cond*/,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
output_mu[i] = viscosity_[regionIdx](p[i]);
|
||||
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
|
||||
}
|
||||
std::fill(output_dmudr, output_dmudr + n, 0.0);
|
||||
}
|
||||
|
||||
void PvtDeadSpline::B(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_B) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
output_B[i] = 1.0/b_[regionIdx](p[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void PvtDeadSpline::dBdp(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_B,
|
||||
double* output_dBdp) const
|
||||
{
|
||||
B(n, p, 0, output_B);
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
double Bg = output_B[i];
|
||||
output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void PvtDeadSpline::b(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
output_b[i] = b_[regionIdx](p[i]);
|
||||
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
|
||||
}
|
||||
std::fill(output_dbdr, output_dbdr + n, 0.0);
|
||||
}
|
||||
|
||||
void PvtDeadSpline::b(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
const PhasePresence* /*cond*/,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int regionIdx = pvtTableIdx_[i];
|
||||
output_b[i] = b_[regionIdx](p[i]);
|
||||
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
|
||||
}
|
||||
std::fill(output_dbdr, output_dbdr + n, 0.0);
|
||||
}
|
||||
|
||||
void PvtDeadSpline::rsSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rsSat,
|
||||
double* output_drsSatdp) const
|
||||
{
|
||||
std::fill(output_rsSat, output_rsSat + n, 0.0);
|
||||
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
void PvtDeadSpline::rvSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rvSat,
|
||||
double* output_drvSatdp) const
|
||||
{
|
||||
std::fill(output_rvSat, output_rvSat + n, 0.0);
|
||||
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
void PvtDeadSpline::R(const int n,
|
||||
const double* /*p*/,
|
||||
const double* /*z*/,
|
||||
double* output_R) const
|
||||
{
|
||||
std::fill(output_R, output_R + n, 0.0);
|
||||
}
|
||||
|
||||
void PvtDeadSpline::dRdp(const int n,
|
||||
const double* /*p*/,
|
||||
const double* /*z*/,
|
||||
double* output_R,
|
||||
double* output_dRdp) const
|
||||
{
|
||||
std::fill(output_R, output_R + n, 0.0);
|
||||
std::fill(output_dRdp, output_dRdp + n, 0.0);
|
||||
}
|
||||
|
||||
}
|
||||
@ -17,10 +17,10 @@
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
|
||||
#define OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
|
||||
#ifndef OPM_PVTDEADSPLINE_HEADER_INCLUDED
|
||||
#define OPM_PVTDEADSPLINE_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
|
||||
#include <opm/core/props/pvt/PvtInterface.hpp>
|
||||
#include <opm/core/utility/UniformTableLinear.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
|
||||
@ -38,12 +38,19 @@ namespace Opm
|
||||
/// are expected to be of size n, size n and n*num_phases, respectively.
|
||||
/// Output arrays shall be of size n, and must be valid before
|
||||
/// calling the method.
|
||||
class SinglePvtDeadSpline : public SinglePvtInterface
|
||||
class PvtDeadSpline : public PvtInterface
|
||||
{
|
||||
public:
|
||||
SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples);
|
||||
SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples);
|
||||
virtual ~SinglePvtDeadSpline();
|
||||
PvtDeadSpline();
|
||||
|
||||
void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
|
||||
const std::vector<int> &pvtTableIdx,
|
||||
int numSamples);
|
||||
void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
|
||||
const std::vector<int> &pvtTableIdx,
|
||||
int numSamples);
|
||||
|
||||
virtual ~PvtDeadSpline();
|
||||
|
||||
/// Viscosity as a function of p and z.
|
||||
virtual void mu(const int n,
|
||||
@ -128,11 +135,12 @@ namespace Opm
|
||||
double* output_dRdp) const;
|
||||
private:
|
||||
// PVT properties of dry gas or dead oil
|
||||
UniformTableLinear<double> b_;
|
||||
UniformTableLinear<double> viscosity_;
|
||||
std::vector<int> pvtTableIdx_;
|
||||
std::vector<UniformTableLinear<double> > b_;
|
||||
std::vector<UniformTableLinear<double> > viscosity_;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
|
||||
#endif // OPM_PVTDEADSPLINE_HEADER_INCLUDED
|
||||
|
||||
79
opm/core/props/pvt/PvtInterface.cpp
Normal file
79
opm/core/props/pvt/PvtInterface.cpp
Normal file
@ -0,0 +1,79 @@
|
||||
/*
|
||||
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/props/pvt/PvtInterface.hpp>
|
||||
|
||||
#include <cassert>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
PvtInterface::PvtInterface()
|
||||
: num_phases_(MaxNumPhases)
|
||||
{
|
||||
for (int i = 0; i < MaxNumPhases; ++i) {
|
||||
phase_pos_[i] = i;
|
||||
}
|
||||
}
|
||||
|
||||
PvtInterface::~PvtInterface()
|
||||
{
|
||||
}
|
||||
|
||||
void PvtInterface::setPhaseConfiguration(const int num_phases, const int* phase_pos)
|
||||
{
|
||||
num_phases_ = num_phases;
|
||||
for (int i = 0; i < MaxNumPhases; ++i) {
|
||||
phase_pos_[i] = phase_pos[i];
|
||||
}
|
||||
}
|
||||
|
||||
void extractPvtTableIndex(std::vector<int> &pvtTableIdx,
|
||||
Opm::DeckConstPtr deck,
|
||||
size_t numCompressed,
|
||||
const int *compressedToCartesianCellIdx)
|
||||
{
|
||||
if (deck->hasKeyword("PVTNUM")) {
|
||||
// we have a "real" multi-pvt deck.
|
||||
|
||||
// get the PVT number from the deck
|
||||
Opm::DeckKeywordConstPtr pvtnumKeyword = deck->getKeyword("PVTNUM");
|
||||
const std::vector<int> &pvtnumData = pvtnumKeyword->getIntData();
|
||||
|
||||
// convert this into an array of compressed cells. since
|
||||
// Eclipse uses Fortran-style indices which start at 1
|
||||
// instead of 0, we subtract 1.
|
||||
pvtTableIdx.resize(numCompressed);
|
||||
for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) {
|
||||
size_t cartesianCellIdx = compressedToCartesianCellIdx[cellIdx];
|
||||
assert(cartesianCellIdx < pvtnumData.size());
|
||||
|
||||
pvtTableIdx[cellIdx] = pvtnumData[cartesianCellIdx] - 1;
|
||||
}
|
||||
}
|
||||
else {
|
||||
// create the pvtIdxArray: all cells use the one and only
|
||||
// PVT table...
|
||||
pvtTableIdx.resize(numCompressed);
|
||||
std::fill(pvtTableIdx.begin(), pvtTableIdx.end(), 0);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
@ -17,22 +17,22 @@
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED
|
||||
#define OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED
|
||||
|
||||
#ifndef OPM_PVTINTERFACE_HEADER_INCLUDED
|
||||
#define OPM_PVTINTERFACE_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
class SinglePvtInterface : public BlackoilPhases
|
||||
class PvtInterface : public BlackoilPhases
|
||||
{
|
||||
public:
|
||||
SinglePvtInterface();
|
||||
PvtInterface();
|
||||
|
||||
virtual ~SinglePvtInterface();
|
||||
virtual ~PvtInterface();
|
||||
|
||||
/// \param[in] num_phases The number of active phases.
|
||||
/// \param[in] phase_pos Array of BlackpoilPhases::MaxNumPhases
|
||||
@ -138,7 +138,22 @@ namespace Opm
|
||||
int phase_pos_[MaxNumPhases];
|
||||
};
|
||||
|
||||
/*!
|
||||
* \brief Helper function to create an array containing the (C-Style)
|
||||
* PVT table index for each compressed cell.
|
||||
*
|
||||
* The main point of this function is to avoid code duplication
|
||||
* because the Eclipse deck only contains Fortran-style PVT table
|
||||
* indices which start at 1 instead of 0 and -- more relevantly -- it
|
||||
* uses logically cartesian cell indices to specify the table index of
|
||||
* a cell.
|
||||
*/
|
||||
void extractPvtTableIndex(std::vector<int> &pvtTableIdx,
|
||||
Opm::DeckConstPtr deck,
|
||||
size_t numCompressed,
|
||||
const int *compressedToCartesianIdx);
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED
|
||||
#endif // OPM_PVTINTERFACE_HEADER_INCLUDED
|
||||
|
||||
516
opm/core/props/pvt/PvtLiveGas.cpp
Normal file
516
opm/core/props/pvt/PvtLiveGas.cpp
Normal file
@ -0,0 +1,516 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: MiscibilityLiveGas.cpp
|
||||
//
|
||||
// Created: Wed Feb 10 09:21:53 2010
|
||||
//
|
||||
// Author: Bjørn Spjelkavik <bsp@sintef.no>
|
||||
//
|
||||
// Revision: $Id$
|
||||
//
|
||||
//===========================================================================
|
||||
/*
|
||||
Copyright 2010 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/pvt/PvtLiveGas.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/linearInterpolation.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
using Opm::linearInterpolation;
|
||||
using Opm::linearInterpolationDerivative;
|
||||
|
||||
|
||||
//------------------------------------------------------------------------
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
PvtLiveGas::PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword,
|
||||
const std::vector<int> &pvtTableIdx)
|
||||
{
|
||||
pvtTableIdx_ = pvtTableIdx;
|
||||
|
||||
int numTables = Opm::PvtgTable::numTables(pvtgKeyword);
|
||||
saturated_gas_table_.resize(numTables);
|
||||
undersat_gas_tables_.resize(numTables);
|
||||
|
||||
for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) {
|
||||
Opm::PvtgTable pvtgTable(pvtgKeyword, pvtTableIdx);
|
||||
|
||||
// GAS, PVTG
|
||||
saturated_gas_table_[pvtTableIdx].resize(4);
|
||||
saturated_gas_table_[pvtTableIdx][0] = pvtgTable.getOuterTable()->getPressureColumn();
|
||||
saturated_gas_table_[pvtTableIdx][1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn();
|
||||
saturated_gas_table_[pvtTableIdx][2] = pvtgTable.getOuterTable()->getGasViscosityColumn();
|
||||
saturated_gas_table_[pvtTableIdx][3] = pvtgTable.getOuterTable()->getOilSolubilityColumn();
|
||||
|
||||
int sz = pvtgTable.getOuterTable()->numRows();
|
||||
undersat_gas_tables_[pvtTableIdx].resize(sz);
|
||||
for (int i=0; i<sz; ++i) {
|
||||
const auto &undersatTable = *pvtgTable.getInnerTable(i);
|
||||
|
||||
undersat_gas_tables_[pvtTableIdx][i].resize(3);
|
||||
undersat_gas_tables_[pvtTableIdx][i][0] = undersatTable.getOilSolubilityColumn();
|
||||
undersat_gas_tables_[pvtTableIdx][i][1] = undersatTable.getGasFormationFactorColumn();
|
||||
undersat_gas_tables_[pvtTableIdx][i][2] = undersatTable.getGasViscosityColumn();
|
||||
}
|
||||
|
||||
// Bg -> 1/Bg
|
||||
for (int i=0; i<sz; ++i) {
|
||||
saturated_gas_table_[pvtTableIdx][1][i] = 1.0/saturated_gas_table_[pvtTableIdx][1][i];
|
||||
for (size_t j=0; j<undersat_gas_tables_[pvtTableIdx][i][1].size(); ++j) {
|
||||
undersat_gas_tables_[pvtTableIdx][i][1][j] = 1.0/undersat_gas_tables_[pvtTableIdx][i][1][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Destructor
|
||||
PvtLiveGas::~PvtLiveGas()
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
void PvtLiveGas::mu(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_mu) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_mu[i] = miscible_gas(p[i], z + num_phases_*i, pvtTableIdx_[i], 2, false);
|
||||
}
|
||||
}
|
||||
|
||||
/// Viscosity and its derivatives as a function of p and r.
|
||||
void PvtLiveGas::mu(const int /*n*/,
|
||||
const double* /*p*/,
|
||||
const double* /*r*/,
|
||||
double* /*output_mu*/,
|
||||
double* /*output_dmudp*/,
|
||||
double* /*output_dmudr*/) const
|
||||
{
|
||||
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
|
||||
}
|
||||
|
||||
/// Viscosity and its derivatives as a function of p and r.
|
||||
void PvtLiveGas::mu(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
const PhasePresence* cond,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const PhasePresence& cnd = cond[i];
|
||||
output_mu[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 2, 0);
|
||||
output_dmudp[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 2, 1);
|
||||
output_dmudr[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 2, 2);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/// Formation volume factor as a function of p and z.
|
||||
void PvtLiveGas::B(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_B) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_B[i] = evalB(p[i], z + num_phases_*i, pvtTableIdx_[i]);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/// Formation volume factor and p-derivative as functions of p and z.
|
||||
void PvtLiveGas::dBdp(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_B,
|
||||
double* output_dBdp) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
evalBDeriv(p[i], z + num_phases_*i, pvtTableIdx_[i], output_B[i], output_dBdp[i]);
|
||||
}
|
||||
}
|
||||
|
||||
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
|
||||
void PvtLiveGas::b(const int /*n*/,
|
||||
const double* /*p*/,
|
||||
const double* /*r*/,
|
||||
double* /*output_b*/,
|
||||
double* /*output_dbdp*/,
|
||||
double* /*output_dbdr*/) const
|
||||
|
||||
{
|
||||
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
|
||||
}
|
||||
|
||||
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
|
||||
void PvtLiveGas::b(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
const PhasePresence* cond,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const PhasePresence& cnd = cond[i];
|
||||
|
||||
output_b[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 1, 0);
|
||||
output_dbdp[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 1, 1);
|
||||
output_dbdr[i] = miscible_gas(p[i], r[i], cnd, pvtTableIdx_[i], 1, 2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
/// Gas resolution and its derivatives at bublepoint as a function of p.
|
||||
void PvtLiveGas::rvSat(const int n,
|
||||
const double* p,
|
||||
double* output_rvSat,
|
||||
double* output_drvSatdp) const
|
||||
{
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
output_rvSat[i] = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
|
||||
saturated_gas_table_[pvtTableIdx][3],p[i]);
|
||||
output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0],
|
||||
saturated_gas_table_[pvtTableIdx][3],p[i]);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
void PvtLiveGas::rsSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rsSat,
|
||||
double* output_drsSatdp) const
|
||||
{
|
||||
std::fill(output_rsSat, output_rsSat + n, 0.0);
|
||||
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
/// Solution factor as a function of p and z.
|
||||
void PvtLiveGas::R(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_R) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_R[i] = evalR(p[i], z + num_phases_*i, pvtTableIdx_[i]);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/// Solution factor and p-derivative as functions of p and z.
|
||||
void PvtLiveGas::dRdp(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_R,
|
||||
double* output_dRdp) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
evalRDeriv(p[i], z + num_phases_*i, pvtTableIdx_[i], output_R[i], output_dRdp[i]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ---- Private methods ----
|
||||
|
||||
double PvtLiveGas::evalB(const double press, const double* surfvol, int pvtTableIdx) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Vapour]] == 0.0) {
|
||||
// To handle no-gas case.
|
||||
return 1.0;
|
||||
}
|
||||
return 1.0/miscible_gas(press, surfvol, pvtTableIdx, 1, false);
|
||||
}
|
||||
|
||||
void PvtLiveGas::evalBDeriv(const double press, const double* surfvol, int pvtTableIdx,
|
||||
double& Bval, double& dBdpval) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Vapour]] == 0.0) {
|
||||
// To handle no-gas case.
|
||||
Bval = 1.0;
|
||||
dBdpval = 0.0;
|
||||
return;
|
||||
}
|
||||
Bval = evalB(press, surfvol, pvtTableIdx);
|
||||
dBdpval = -Bval*Bval*miscible_gas(press, surfvol, pvtTableIdx, 1, true);
|
||||
}
|
||||
|
||||
double PvtLiveGas::evalR(const double press, const double* surfvol, int pvtTableIdx) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Liquid]] == 0.0) {
|
||||
// To handle no-gas case.
|
||||
return 0.0;
|
||||
}
|
||||
double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
|
||||
saturated_gas_table_[pvtTableIdx][3], press);
|
||||
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
|
||||
if (satR < maxR ) {
|
||||
// Saturated case
|
||||
return satR;
|
||||
} else {
|
||||
// Undersaturated case
|
||||
return maxR;
|
||||
}
|
||||
}
|
||||
|
||||
void PvtLiveGas::evalRDeriv(const double press, const double* surfvol, int pvtTableIdx,
|
||||
double& Rval, double& dRdpval) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Liquid]] == 0.0) {
|
||||
// To handle no-gas case.
|
||||
Rval = 0.0;
|
||||
dRdpval = 0.0;
|
||||
return;
|
||||
}
|
||||
double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
|
||||
saturated_gas_table_[pvtTableIdx][3], press);
|
||||
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
|
||||
if (satR < maxR ) {
|
||||
// Saturated case
|
||||
Rval = satR;
|
||||
dRdpval = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0],
|
||||
saturated_gas_table_[pvtTableIdx][3],
|
||||
press);
|
||||
} else {
|
||||
// Undersaturated case
|
||||
Rval = maxR;
|
||||
dRdpval = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
double PvtLiveGas::miscible_gas(const double press,
|
||||
const double* surfvol,
|
||||
const int pvtTableIdx,
|
||||
const int item,
|
||||
const bool deriv) const
|
||||
{
|
||||
const std::vector<std::vector<double> > &saturatedGasTable =
|
||||
saturated_gas_table_[pvtTableIdx];
|
||||
const std::vector<std::vector<std::vector<double> > > &undersatGasTables =
|
||||
undersat_gas_tables_[pvtTableIdx];
|
||||
|
||||
int section;
|
||||
double Rval = linearInterpolation(saturatedGasTable[0],
|
||||
saturatedGasTable[3], press,
|
||||
section);
|
||||
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
|
||||
if (deriv) {
|
||||
if (Rval < maxR ) { // Saturated case
|
||||
return linearInterpolationDerivative(saturatedGasTable[0],
|
||||
saturatedGasTable[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = section;
|
||||
if (undersatGasTables[is][0].size() < 2) {
|
||||
double val = (saturatedGasTable[item][is+1]
|
||||
- saturatedGasTable[item][is]) /
|
||||
(saturatedGasTable[0][is+1] -
|
||||
saturatedGasTable[0][is]);
|
||||
return val;
|
||||
}
|
||||
double val1 =
|
||||
linearInterpolation(undersatGasTables[is][0],
|
||||
undersatGasTables[is][item],
|
||||
maxR);
|
||||
double val2 =
|
||||
linearInterpolation(undersatGasTables[is+1][0],
|
||||
undersatGasTables[is+1][item],
|
||||
maxR);
|
||||
double val = (val2 - val1)/
|
||||
(saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
|
||||
return val;
|
||||
}
|
||||
} else {
|
||||
if (Rval < maxR ) { // Saturated case
|
||||
return linearInterpolation(saturatedGasTable[0],
|
||||
saturatedGasTable[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = section;
|
||||
// Extrapolate from first table section
|
||||
if (is == 0 && press < saturatedGasTable[0][0]) {
|
||||
return linearInterpolation(undersatGasTables[0][0],
|
||||
undersatGasTables[0][item],
|
||||
maxR);
|
||||
}
|
||||
|
||||
// Extrapolate from last table section
|
||||
int ltp = saturatedGasTable[0].size() - 1;
|
||||
if (is+1 == ltp && press > saturatedGasTable[0][ltp]) {
|
||||
return linearInterpolation(undersatGasTables[ltp][0],
|
||||
undersatGasTables[ltp][item],
|
||||
maxR);
|
||||
}
|
||||
|
||||
// Interpolate between table sections
|
||||
double w = (press - saturatedGasTable[0][is]) /
|
||||
(saturatedGasTable[0][is+1] -
|
||||
saturatedGasTable[0][is]);
|
||||
if (undersatGasTables[is][0].size() < 2) {
|
||||
double val = saturatedGasTable[item][is] +
|
||||
w*(saturatedGasTable[item][is+1] -
|
||||
saturatedGasTable[item][is]);
|
||||
return val;
|
||||
}
|
||||
double val1 =
|
||||
linearInterpolation(undersatGasTables[is][0],
|
||||
undersatGasTables[is][item],
|
||||
maxR);
|
||||
double val2 =
|
||||
linearInterpolation(undersatGasTables[is+1][0],
|
||||
undersatGasTables[is+1][item],
|
||||
maxR);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double PvtLiveGas::miscible_gas(const double press,
|
||||
const double r,
|
||||
const PhasePresence& cond,
|
||||
const int pvtTableIdx,
|
||||
const int item,
|
||||
const int deriv) const
|
||||
{
|
||||
const std::vector<std::vector<double> > &saturatedGasTable =
|
||||
saturated_gas_table_[pvtTableIdx];
|
||||
const std::vector<std::vector<std::vector<double> > > &undersatGasTables =
|
||||
undersat_gas_tables_[pvtTableIdx];
|
||||
|
||||
const bool isSat = cond.hasFreeOil();
|
||||
|
||||
// Derivative w.r.t p
|
||||
if (deriv == 1) {
|
||||
if (isSat) { // Saturated case
|
||||
return linearInterpolationDerivative(saturatedGasTable[0],
|
||||
saturatedGasTable[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturatedGasTable[0], press);
|
||||
if (undersatGasTables[is][0].size() < 2) {
|
||||
double val = (saturatedGasTable[item][is+1]
|
||||
- saturatedGasTable[item][is]) /
|
||||
(saturatedGasTable[0][is+1] -
|
||||
saturatedGasTable[0][is]);
|
||||
return val;
|
||||
}
|
||||
double val1 =
|
||||
linearInterpolation(undersatGasTables[is][0],
|
||||
undersatGasTables[is][item],
|
||||
r);
|
||||
double val2 =
|
||||
linearInterpolation(undersatGasTables[is+1][0],
|
||||
undersatGasTables[is+1][item],
|
||||
r);
|
||||
double val = (val2 - val1)/
|
||||
(saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
|
||||
return val;
|
||||
}
|
||||
} else if (deriv == 2){
|
||||
if (isSat) {
|
||||
return 0;
|
||||
} else {
|
||||
int is = tableIndex(saturatedGasTable[0], press);
|
||||
double w = (press - saturatedGasTable[0][is]) /
|
||||
(saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
|
||||
assert(undersatGasTables[is][0].size() >= 2);
|
||||
assert(undersatGasTables[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolationDerivative(undersatGasTables[is][0],
|
||||
undersatGasTables[is][item],
|
||||
r);
|
||||
double val2 =
|
||||
linearInterpolationDerivative(undersatGasTables[is+1][0],
|
||||
undersatGasTables[is+1][item],
|
||||
r);
|
||||
|
||||
double val = val1 + w * (val2 - val1);
|
||||
return val;
|
||||
|
||||
}
|
||||
} else {
|
||||
if (isSat) { // Saturated case
|
||||
return linearInterpolation(saturatedGasTable[0],
|
||||
saturatedGasTable[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturatedGasTable[0], press);
|
||||
// Extrapolate from first table section
|
||||
if (is == 0 && press < saturatedGasTable[0][0]) {
|
||||
return linearInterpolation(undersatGasTables[0][0],
|
||||
undersatGasTables[0][item],
|
||||
r);
|
||||
}
|
||||
|
||||
// Extrapolate from last table section
|
||||
//int ltp = saturatedGasTable[0].size() - 1;
|
||||
//if (is+1 == ltp && press > saturatedGasTable[0][ltp]) {
|
||||
// return linearInterpolation(undersatGasTables[ltp][0],
|
||||
// undersatGasTables[ltp][item],
|
||||
// r);
|
||||
//}
|
||||
|
||||
// Interpolate between table sections
|
||||
double w = (press - saturatedGasTable[0][is]) /
|
||||
(saturatedGasTable[0][is+1] -
|
||||
saturatedGasTable[0][is]);
|
||||
if (undersatGasTables[is][0].size() < 2) {
|
||||
double val = saturatedGasTable[item][is] +
|
||||
w*(saturatedGasTable[item][is+1] -
|
||||
saturatedGasTable[item][is]);
|
||||
return val;
|
||||
}
|
||||
double val1 =
|
||||
linearInterpolation(undersatGasTables[is][0],
|
||||
undersatGasTables[is][item],
|
||||
r);
|
||||
double val2 =
|
||||
linearInterpolation(undersatGasTables[is+1][0],
|
||||
undersatGasTables[is+1][item],
|
||||
r);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
@ -17,10 +17,10 @@
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
|
||||
#define OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
|
||||
#ifndef OPM_PVTLIVEGAS_HEADER_INCLUDED
|
||||
#define OPM_PVTLIVEGAS_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
|
||||
#include <opm/core/props/pvt/PvtInterface.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/PvtgTable.hpp>
|
||||
|
||||
@ -35,11 +35,11 @@ namespace Opm
|
||||
/// are expected to be of size n, size n and n*num_phases, respectively.
|
||||
/// Output arrays shall be of size n, and must be valid before
|
||||
/// calling the method.
|
||||
class SinglePvtLiveGas : public SinglePvtInterface
|
||||
class PvtLiveGas : public PvtInterface
|
||||
{
|
||||
public:
|
||||
SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable);
|
||||
virtual ~SinglePvtLiveGas();
|
||||
PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword, const std::vector<int> &pvtTableIdx);
|
||||
virtual ~PvtLiveGas();
|
||||
|
||||
/// Viscosity as a function of p and z.
|
||||
virtual void mu(const int n,
|
||||
@ -126,28 +126,30 @@ namespace Opm
|
||||
double* output_dRdp) const;
|
||||
|
||||
protected:
|
||||
double evalB(double press, const double* surfvol) const;
|
||||
void evalBDeriv(double press, const double* surfvol, double& B, double& dBdp) const;
|
||||
double evalR(double press, const double* surfvol) const;
|
||||
void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const;
|
||||
double evalB(double press, const double* surfvol, int pvtTableIdx) const;
|
||||
void evalBDeriv(double press, const double* surfvol, int pvtTableIdx, double& B, double& dBdp) const;
|
||||
double evalR(double press, const double* surfvol, int pvtTableIdx) const;
|
||||
void evalRDeriv(double press, const double* surfvol, int pvtTableIdx, double& R, double& dRdp) const;
|
||||
|
||||
// item: 1=>B 2=>mu;
|
||||
double miscible_gas(const double press,
|
||||
const double* surfvol,
|
||||
const int pvtTableIdx,
|
||||
const int item,
|
||||
const bool deriv = false) const;
|
||||
double miscible_gas(const double press,
|
||||
const double r,
|
||||
const PhasePresence& cond,
|
||||
const int pvtTableIdx,
|
||||
const int item,
|
||||
const int deriv = 0) const;
|
||||
// PVT properties of wet gas (with vaporised oil)
|
||||
std::vector<std::vector<double> > saturated_gas_table_;
|
||||
std::vector<std::vector<std::vector<double> > > undersat_gas_tables_;
|
||||
|
||||
std::vector<int> pvtTableIdx_;
|
||||
std::vector< std::vector<std::vector<double> > > saturated_gas_table_;
|
||||
std::vector< std::vector<std::vector<std::vector<double> > > > undersat_gas_tables_;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
|
||||
#endif // OPM_PVTLIVEGAS_HEADER_INCLUDED
|
||||
|
||||
583
opm/core/props/pvt/PvtLiveOil.cpp
Normal file
583
opm/core/props/pvt/PvtLiveOil.cpp
Normal file
@ -0,0 +1,583 @@
|
||||
/*
|
||||
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/props/pvt/PvtLiveOil.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/linearInterpolation.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
using Opm::linearInterpolation;
|
||||
using Opm::linearInterpolationDerivative;
|
||||
using Opm::tableIndex;
|
||||
|
||||
//------------------------------------------------------------------------
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
PvtLiveOil::PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword,
|
||||
const std::vector<int> &pvtTableIdx)
|
||||
{
|
||||
pvtTableIdx_ = pvtTableIdx;
|
||||
|
||||
int numTables = Opm::PvtoTable::numTables(pvtoKeyword);
|
||||
saturated_oil_table_.resize(numTables);
|
||||
undersat_oil_tables_.resize(numTables);
|
||||
|
||||
for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) {
|
||||
Opm::PvtoTable pvtoTable(pvtoKeyword, pvtTableIdx);
|
||||
|
||||
const auto saturatedPvto = pvtoTable.getOuterTable();
|
||||
|
||||
// OIL, PVTO
|
||||
saturated_oil_table_[pvtTableIdx].resize(4);
|
||||
const int sz = saturatedPvto->numRows();
|
||||
for (int k=0; k<4; ++k) {
|
||||
saturated_oil_table_[pvtTableIdx][k].resize(sz);
|
||||
}
|
||||
for (int i=0; i<sz; ++i) {
|
||||
saturated_oil_table_[pvtTableIdx][0][i] = saturatedPvto->getPressureColumn()[i]; // p
|
||||
saturated_oil_table_[pvtTableIdx][1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo
|
||||
saturated_oil_table_[pvtTableIdx][2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o
|
||||
saturated_oil_table_[pvtTableIdx][3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs
|
||||
}
|
||||
|
||||
undersat_oil_tables_[pvtTableIdx].resize(sz);
|
||||
for (int i=0; i<sz; ++i) {
|
||||
const auto undersaturatedPvto = pvtoTable.getInnerTable(i);
|
||||
|
||||
undersat_oil_tables_[pvtTableIdx][i].resize(3);
|
||||
int tsize = undersaturatedPvto->numRows();
|
||||
undersat_oil_tables_[pvtTableIdx][i][0].resize(tsize);
|
||||
undersat_oil_tables_[pvtTableIdx][i][1].resize(tsize);
|
||||
undersat_oil_tables_[pvtTableIdx][i][2].resize(tsize);
|
||||
for (int j=0; j<tsize; ++j) {
|
||||
undersat_oil_tables_[pvtTableIdx][i][0][j] = undersaturatedPvto->getPressureColumn()[j]; // p
|
||||
undersat_oil_tables_[pvtTableIdx][i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j]; // 1/Bo
|
||||
undersat_oil_tables_[pvtTableIdx][i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j]; // mu_o
|
||||
}
|
||||
}
|
||||
|
||||
// Complete undersaturated tables by extrapolating from existing data
|
||||
// as is done in Eclipse and Mrst
|
||||
int iNext = -1;
|
||||
for (int i=0; i<sz; ++i) {
|
||||
// Skip records already containing undersaturated data
|
||||
if (undersat_oil_tables_[pvtTableIdx][i][0].size() > 1) {
|
||||
continue;
|
||||
}
|
||||
// Look ahead for next record containing undersaturated data
|
||||
if (iNext < i) {
|
||||
iNext = i+1;
|
||||
while (iNext<sz && undersat_oil_tables_[pvtTableIdx][iNext][0].size() < 2) {
|
||||
++iNext;
|
||||
}
|
||||
if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
|
||||
}
|
||||
// Add undersaturated data to current record while maintaining compressibility and viscosibility
|
||||
typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
|
||||
for (sz_t j=1; j<undersat_oil_tables_[pvtTableIdx][iNext][0].size(); ++j) {
|
||||
double diffPressure = undersat_oil_tables_[pvtTableIdx][iNext][0][j]-undersat_oil_tables_[pvtTableIdx][iNext][0][j-1];
|
||||
double pressure = undersat_oil_tables_[pvtTableIdx][i][0].back()+diffPressure;
|
||||
undersat_oil_tables_[pvtTableIdx][i][0].push_back(pressure);
|
||||
double compr = (1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j]-1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j-1])
|
||||
/ (0.5*(1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j]+1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j-1]));
|
||||
double B = (1.0/undersat_oil_tables_[pvtTableIdx][i][1].back())*(1.0+0.5*compr)/(1.0-0.5*compr);
|
||||
undersat_oil_tables_[pvtTableIdx][i][1].push_back(1.0/B);
|
||||
double visc = (undersat_oil_tables_[pvtTableIdx][iNext][2][j]-undersat_oil_tables_[pvtTableIdx][iNext][2][j-1])
|
||||
/ (0.5*(undersat_oil_tables_[pvtTableIdx][iNext][2][j]+undersat_oil_tables_[pvtTableIdx][iNext][2][j-1]));
|
||||
double mu = (undersat_oil_tables_[pvtTableIdx][i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc);
|
||||
undersat_oil_tables_[pvtTableIdx][i][2].push_back(mu);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Destructor.
|
||||
PvtLiveOil::~PvtLiveOil()
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
/// Viscosity as a function of p and z.
|
||||
void PvtLiveOil::mu(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_mu) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
|
||||
output_mu[i] = miscible_oil(p[i], z + num_phases_*i, pvtTableIdx, 2, false);
|
||||
}
|
||||
}
|
||||
|
||||
/// Viscosity and its derivatives as a function of p and r.
|
||||
void PvtLiveOil::mu(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
|
||||
output_mu[i] = miscible_oil(p[i], r[i], pvtTableIdx, 2, 0);
|
||||
output_dmudp[i] = miscible_oil(p[i], r[i], pvtTableIdx, 2, 1);
|
||||
output_dmudr[i] = miscible_oil(p[i], r[i], pvtTableIdx, 2, 2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
/// Viscosity and its derivatives as a function of p and r.
|
||||
void PvtLiveOil::mu(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
const PhasePresence* cond,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
const PhasePresence& cnd = cond[i];
|
||||
|
||||
output_mu[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 2, 0);
|
||||
output_dmudp[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 2, 1);
|
||||
output_dmudr[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 2, 2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// Formation volume factor as a function of p and z.
|
||||
void PvtLiveOil::B(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_B) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
|
||||
output_B[i] = evalB(pvtTableIdx, p[i], z + num_phases_*i);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/// Formation volume factor and p-derivative as functions of p and z.
|
||||
void PvtLiveOil::dBdp(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_B,
|
||||
double* output_dBdp) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
|
||||
evalBDeriv(pvtTableIdx, p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void PvtLiveOil::b(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
|
||||
output_b[i] = miscible_oil(pvtTableIdx, p[i], r[i], 1, 0);
|
||||
output_dbdp[i] = miscible_oil(pvtTableIdx, p[i], r[i], 1, 1);
|
||||
output_dbdr[i] = miscible_oil(pvtTableIdx, p[i], r[i], 1, 2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
void PvtLiveOil::b(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
const PhasePresence* cond,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const PhasePresence& cnd = cond[i];
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
|
||||
output_b[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 1, 0);
|
||||
output_dbdp[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 1, 1);
|
||||
output_dbdr[i] = miscible_oil(p[i], r[i], cnd, pvtTableIdx, 1, 2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
void PvtLiveOil::rsSat(const int n,
|
||||
const double* p,
|
||||
double* output_rsSat,
|
||||
double* output_drsSatdp) const
|
||||
{
|
||||
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
output_rsSat[i] = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][3],p[i]);
|
||||
output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][3],p[i]);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
void PvtLiveOil::rvSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rvSat,
|
||||
double* output_drvSatdp) const
|
||||
{
|
||||
std::fill(output_rvSat, output_rvSat + n, 0.0);
|
||||
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
/// Solution factor as a function of p and z.
|
||||
void PvtLiveOil::R(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_R) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
output_R[i] = evalR(pvtTableIdx, p[i], z + num_phases_*i);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/// Solution factor and p-derivative as functions of p and z.
|
||||
void PvtLiveOil::dRdp(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_R,
|
||||
double* output_dRdp) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
int pvtTableIdx = pvtTableIdx_[i];
|
||||
evalRDeriv(pvtTableIdx, p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
// ---- Private methods ----
|
||||
|
||||
double PvtLiveOil::evalB(size_t pvtTableIdx,
|
||||
double press, const double* surfvol) const
|
||||
{
|
||||
// if (surfvol[phase_pos_[Liquid]] == 0.0) return 1.0; // To handle no-oil case.
|
||||
return 1.0/miscible_oil(press, surfvol, pvtTableIdx, 1, false);
|
||||
}
|
||||
|
||||
|
||||
void PvtLiveOil::evalBDeriv(size_t pvtTableIdx,
|
||||
const double press, const double* surfvol,
|
||||
double& Bval, double& dBdpval) const
|
||||
{
|
||||
Bval = evalB(pvtTableIdx, press, surfvol);
|
||||
dBdpval = -Bval*Bval*miscible_oil(press, surfvol, pvtTableIdx, 1, true);
|
||||
}
|
||||
|
||||
double PvtLiveOil::evalR(size_t pvtTableIdx,
|
||||
double press, const double* surfvol) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Vapour]] == 0.0) {
|
||||
return 0.0;
|
||||
}
|
||||
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][3], press);
|
||||
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
|
||||
if (Rval < maxR ) { // Saturated case
|
||||
return Rval;
|
||||
} else {
|
||||
return maxR; // Undersaturated case
|
||||
}
|
||||
}
|
||||
|
||||
void PvtLiveOil::evalRDeriv(size_t pvtTableIdx,
|
||||
const double press, const double* surfvol,
|
||||
double& Rval, double& dRdpval) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Vapour]] == 0.0) {
|
||||
Rval = 0.0;
|
||||
dRdpval = 0.0;
|
||||
return;
|
||||
}
|
||||
Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][3], press);
|
||||
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
|
||||
if (Rval < maxR ) {
|
||||
// Saturated case
|
||||
dRdpval = linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][3],
|
||||
press);
|
||||
} else {
|
||||
// Undersaturated case
|
||||
Rval = maxR;
|
||||
dRdpval = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double PvtLiveOil::miscible_oil(const double press,
|
||||
const double* surfvol,
|
||||
const int pvtTableIdx,
|
||||
const int item,
|
||||
const bool deriv) const
|
||||
{
|
||||
int section;
|
||||
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][3],
|
||||
press, section);
|
||||
double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
|
||||
if (deriv) {
|
||||
if (Rval < maxR ) { // Saturated case
|
||||
return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], maxR);
|
||||
double w = (maxR - saturated_oil_table_[pvtTableIdx][3][is]) /
|
||||
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
} else {
|
||||
if (Rval < maxR ) { // Saturated case
|
||||
return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
// Interpolate between table sections
|
||||
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], maxR);
|
||||
double w = (maxR - saturated_oil_table_[pvtTableIdx][3][is]) /
|
||||
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double PvtLiveOil::miscible_oil(const double press,
|
||||
const double r,
|
||||
const int pvtTableIdx,
|
||||
const int item,
|
||||
const int deriv) const
|
||||
{
|
||||
int section;
|
||||
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][3],
|
||||
press, section);
|
||||
// derivative with respect to frist component (pressure)
|
||||
if (deriv == 1) {
|
||||
if (Rval <= r ) { // Saturated case
|
||||
return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
|
||||
double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
|
||||
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
// derivative with respect to second component (r)
|
||||
} else if (deriv == 2) {
|
||||
if (Rval <= r ) { // Saturated case
|
||||
return 0;
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is+1][item],
|
||||
press);
|
||||
|
||||
double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][3][is+1]-saturated_oil_table_[pvtTableIdx][3][is]);
|
||||
return val;
|
||||
}
|
||||
|
||||
|
||||
} else {
|
||||
if (Rval <= r ) { // Saturated case
|
||||
return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
// Interpolate between table sections
|
||||
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
|
||||
double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
|
||||
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double PvtLiveOil::miscible_oil(const double press,
|
||||
const double r,
|
||||
const PhasePresence& cond,
|
||||
const int pvtTableIdx,
|
||||
const int item,
|
||||
const int deriv) const
|
||||
{
|
||||
const bool isSat = cond.hasFreeGas();
|
||||
|
||||
// derivative with respect to frist component (pressure)
|
||||
if (deriv == 1) {
|
||||
if (isSat) { // Saturated case
|
||||
return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
|
||||
double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
|
||||
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
// derivative with respect to second component (r)
|
||||
} else if (deriv == 2) {
|
||||
if (isSat) { // Saturated case
|
||||
return 0;
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is+1][item],
|
||||
press);
|
||||
|
||||
double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][3][is+1]-saturated_oil_table_[pvtTableIdx][3][is]);
|
||||
return val;
|
||||
}
|
||||
|
||||
|
||||
} else {
|
||||
if (isSat) { // Saturated case
|
||||
return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
|
||||
saturated_oil_table_[pvtTableIdx][item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
// Interpolate between table sections
|
||||
int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
|
||||
double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
|
||||
(saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
|
||||
undersat_oil_tables_[pvtTableIdx][is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
@ -17,11 +17,12 @@
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
|
||||
#define OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
|
||||
#ifndef OPM_PVTLIVEOIL_HEADER_INCLUDED
|
||||
#define OPM_PVTLIVEOIL_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
|
||||
#include <opm/core/props/pvt/PvtInterface.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
|
||||
|
||||
#include <vector>
|
||||
@ -35,11 +36,11 @@ namespace Opm
|
||||
/// are expected to be of size n, size n and n*num_phases, respectively.
|
||||
/// Output arrays shall be of size n, and must be valid before
|
||||
/// calling the method.
|
||||
class SinglePvtLiveOil : public SinglePvtInterface
|
||||
class PvtLiveOil : public PvtInterface
|
||||
{
|
||||
public:
|
||||
SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable);
|
||||
virtual ~SinglePvtLiveOil();
|
||||
PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword, const std::vector<int> &pvtTableIdx);
|
||||
virtual ~PvtLiveOil();
|
||||
|
||||
/// Viscosity as a function of p and z.
|
||||
virtual void mu(const int n,
|
||||
@ -124,34 +125,38 @@ namespace Opm
|
||||
double* output_dRdp) const;
|
||||
|
||||
private:
|
||||
double evalB(double press, const double* surfvol) const;
|
||||
void evalBDeriv(double press, const double* surfvol, double& B, double& dBdp) const;
|
||||
double evalR(double press, const double* surfvol) const;
|
||||
void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const;
|
||||
double evalB(size_t pvtTableIdx, double press, const double* surfvol) const;
|
||||
void evalBDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& B, double& dBdp) const;
|
||||
double evalR(size_t pvtTableIdx, double press, const double* surfvol) const;
|
||||
void evalRDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& R, double& dRdp) const;
|
||||
|
||||
// item: 1=>1/B 2=>mu;
|
||||
double miscible_oil(const double press,
|
||||
const double* surfvol,
|
||||
const int pvtTableIdx,
|
||||
const int item,
|
||||
const bool deriv = false) const;
|
||||
|
||||
double miscible_oil(const double press,
|
||||
const double r,
|
||||
const int pvtTableIdx,
|
||||
const int item,
|
||||
const int deriv = 0) const;
|
||||
|
||||
double miscible_oil(const double press,
|
||||
const double r,
|
||||
const PhasePresence& cond,
|
||||
const int pvtTableIdx,
|
||||
const int item,
|
||||
const int deriv = 0) const;
|
||||
|
||||
// PVT properties of live oil (with dissolved gas)
|
||||
std::vector<std::vector<double> > saturated_oil_table_;
|
||||
std::vector<std::vector<std::vector<double> > > undersat_oil_tables_;
|
||||
std::vector<int> pvtTableIdx_;
|
||||
std::vector<std::vector<std::vector<double> > > saturated_oil_table_;
|
||||
std::vector<std::vector<std::vector<std::vector<double> > > > undersat_oil_tables_;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
|
||||
#endif // OPM_PVTLIVEOIL_HEADER_INCLUDED
|
||||
|
||||
@ -35,7 +35,7 @@ namespace Opm
|
||||
|
||||
void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck )
|
||||
{
|
||||
// If we need multiple regions, this class and the SinglePvt* classes must change.
|
||||
// If we need multiple regions, this class and the Pvt* classes must change.
|
||||
int region_number = 0;
|
||||
|
||||
PhaseUsage phase_usage = phaseUsageFromDeck(deck);
|
||||
|
||||
@ -1,291 +0,0 @@
|
||||
/*
|
||||
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
|
||||
#define OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
|
||||
#include <opm/parser/eclipse/Utility/PvcdoTable.hpp>
|
||||
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
/// Class for constant compressible phases (PVTW or PVCDO).
|
||||
/// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
|
||||
/// or pressure (p) and gas resolution factor (r).
|
||||
/// For all the virtual methods, the following apply: p, r and z
|
||||
/// are expected to be of size n, size n and n*num_phases, respectively.
|
||||
/// Output arrays shall be of size n, and must be valid before
|
||||
/// calling the method.
|
||||
class SinglePvtConstCompr : public SinglePvtInterface
|
||||
{
|
||||
public:
|
||||
SinglePvtConstCompr(const Opm::PvtwTable &pvtwTable)
|
||||
{
|
||||
if (pvtwTable.numRows() != 1)
|
||||
OPM_THROW(std::runtime_error,
|
||||
"The table specified by the PVTW keyword is required"
|
||||
"to exhibit exactly one row.");
|
||||
ref_press_ = pvtwTable.getPressureColumn()[0];
|
||||
ref_B_ = pvtwTable.getFormationFactorColumn()[0];
|
||||
comp_ = pvtwTable.getCompressibilityColumn()[0];
|
||||
viscosity_ = pvtwTable.getViscosityColumn()[0];
|
||||
visc_comp_ = pvtwTable.getViscosibilityColumn()[0];
|
||||
}
|
||||
|
||||
SinglePvtConstCompr(const Opm::PvcdoTable &pvcdoTable)
|
||||
{
|
||||
if (pvcdoTable.numRows() != 1)
|
||||
OPM_THROW(std::runtime_error,
|
||||
"The table specified by the PVCDO keyword is required"
|
||||
"to exhibit exactly one row.");
|
||||
ref_press_ = pvcdoTable.getPressureColumn()[0];
|
||||
ref_B_ = pvcdoTable.getFormationFactorColumn()[0];
|
||||
comp_ = pvcdoTable.getCompressibilityColumn()[0];
|
||||
viscosity_ = pvcdoTable.getViscosityColumn()[0];
|
||||
visc_comp_ = pvcdoTable.getViscosibilityColumn()[0];
|
||||
}
|
||||
|
||||
SinglePvtConstCompr(double visc)
|
||||
: ref_press_(0.0),
|
||||
ref_B_(1.0),
|
||||
comp_(0.0),
|
||||
viscosity_(visc),
|
||||
visc_comp_(0.0)
|
||||
{
|
||||
}
|
||||
|
||||
virtual ~SinglePvtConstCompr()
|
||||
{
|
||||
}
|
||||
|
||||
virtual void mu(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_mu) const
|
||||
{
|
||||
if (visc_comp_) {
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
double x = -visc_comp_*(p[i] - ref_press_);
|
||||
output_mu[i] = viscosity_/(1.0 + x + 0.5*x*x);
|
||||
}
|
||||
} else {
|
||||
std::fill(output_mu, output_mu + n, viscosity_);
|
||||
}
|
||||
}
|
||||
|
||||
virtual void mu(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
if (visc_comp_) {
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
double x = -visc_comp_*(p[i] - ref_press_);
|
||||
double d = (1.0 + x + 0.5*x*x);
|
||||
output_mu[i] = viscosity_/d;
|
||||
output_dmudp[i] = (viscosity_/(d*d))*(1+x) * visc_comp_;
|
||||
}
|
||||
} else {
|
||||
std::fill(output_mu, output_mu + n, viscosity_);
|
||||
std::fill(output_dmudp, output_dmudp + n, 0.0);
|
||||
}
|
||||
std::fill(output_dmudr, output_dmudr + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void mu(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
const PhasePresence* /*cond*/,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
if (visc_comp_) {
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
double x = -visc_comp_*(p[i] - ref_press_);
|
||||
double d = (1.0 + x + 0.5*x*x);
|
||||
output_mu[i] = viscosity_/d;
|
||||
output_dmudp[i] = (viscosity_/(d*d))*(1+x) * visc_comp_;
|
||||
}
|
||||
} else {
|
||||
std::fill(output_mu, output_mu + n, viscosity_);
|
||||
std::fill(output_dmudp, output_dmudp + n, 0.0);
|
||||
}
|
||||
std::fill(output_dmudr, output_dmudr + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void B(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_B) const
|
||||
{
|
||||
if (comp_) {
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
double x = comp_*(p[i] - ref_press_);
|
||||
output_B[i] = ref_B_/(1.0 + x + 0.5*x*x);
|
||||
}
|
||||
} else {
|
||||
std::fill(output_B, output_B + n, ref_B_);
|
||||
}
|
||||
}
|
||||
|
||||
virtual void dBdp(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_B,
|
||||
double* output_dBdp) const
|
||||
{
|
||||
if (comp_) {
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
double x = comp_*(p[i] - ref_press_);
|
||||
double d = (1.0 + x + 0.5*x*x);
|
||||
output_B[i] = ref_B_/d;
|
||||
output_dBdp[i] = (-ref_B_/(d*d))*(1 + x) * comp_;
|
||||
}
|
||||
} else {
|
||||
std::fill(output_B, output_B + n, ref_B_);
|
||||
std::fill(output_dBdp, output_dBdp + n, 0.0);
|
||||
}
|
||||
}
|
||||
|
||||
virtual void b(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
{
|
||||
if (comp_) {
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
double x = comp_*(p[i] - ref_press_);
|
||||
double d = (1.0 + x + 0.5*x*x);
|
||||
|
||||
// b = 1/B = d/ref_B_B;
|
||||
output_b[i] = d/ref_B_;
|
||||
output_dbdp[i] = (1 + x) * comp_/ref_B_;
|
||||
}
|
||||
} else {
|
||||
|
||||
std::fill(output_b, output_b + n, 1/ref_B_);
|
||||
std::fill(output_dbdp, output_dbdp + n, 0.0);
|
||||
}
|
||||
std::fill(output_dbdr, output_dbdr + n, 0.0);
|
||||
|
||||
}
|
||||
|
||||
virtual void b(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
const PhasePresence* /*cond*/,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
{
|
||||
if (comp_) {
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
double x = comp_*(p[i] - ref_press_);
|
||||
double d = (1.0 + x + 0.5*x*x);
|
||||
|
||||
// b = 1/B = d/ref_B_B;
|
||||
output_b[i] = d/ref_B_;
|
||||
output_dbdp[i] = (1 + x) * comp_/ref_B_;
|
||||
}
|
||||
} else {
|
||||
|
||||
std::fill(output_b, output_b + n, 1/ref_B_);
|
||||
std::fill(output_dbdp, output_dbdp + n, 0.0);
|
||||
}
|
||||
std::fill(output_dbdr, output_dbdr + n, 0.0);
|
||||
|
||||
}
|
||||
|
||||
virtual void rsSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rsSat,
|
||||
double* output_drsSatdp) const
|
||||
|
||||
{
|
||||
std::fill(output_rsSat, output_rsSat + n, 0.0);
|
||||
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void rvSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rvSat,
|
||||
double* output_drvSatdp) const
|
||||
|
||||
{
|
||||
std::fill(output_rvSat, output_rvSat + n, 0.0);
|
||||
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void R(const int n,
|
||||
const double* /*p*/,
|
||||
const double* /*z*/,
|
||||
double* output_R) const
|
||||
{
|
||||
std::fill(output_R, output_R + n, 0.0);
|
||||
}
|
||||
|
||||
virtual void dRdp(const int n,
|
||||
const double* /*p*/,
|
||||
const double* /*z*/,
|
||||
double* output_R,
|
||||
double* output_dRdp) const
|
||||
{
|
||||
std::fill(output_R, output_R + n, 0.0);
|
||||
std::fill(output_dRdp, output_dRdp + n, 0.0);
|
||||
}
|
||||
|
||||
private:
|
||||
double ref_press_;
|
||||
double ref_B_;
|
||||
double comp_;
|
||||
double viscosity_;
|
||||
double visc_comp_;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif // OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
|
||||
|
||||
@ -1,226 +0,0 @@
|
||||
/*
|
||||
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
|
||||
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/SingleRecordTable.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
// Extra includes for debug dumping of tables.
|
||||
// #include <boost/lexical_cast.hpp>
|
||||
// #include <string>
|
||||
// #include <fstream>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
//------------------------------------------------------------------------
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
|
||||
SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples)
|
||||
{
|
||||
int numRows = pvdoTable.numRows();
|
||||
|
||||
// Copy data
|
||||
const std::vector<double> &press = pvdoTable.getPressureColumn();
|
||||
const std::vector<double> &B = pvdoTable.getFormationFactorColumn();
|
||||
const std::vector<double> &visc = pvdoTable.getViscosityColumn();
|
||||
|
||||
std::vector<double> B_inv(numRows);
|
||||
for (int i = 0; i < numRows; ++i) {
|
||||
B_inv[i] = 1.0 / B[i];
|
||||
}
|
||||
|
||||
buildUniformMonotoneTable(press, B_inv, samples, b_);
|
||||
buildUniformMonotoneTable(press, visc, samples, viscosity_);
|
||||
}
|
||||
|
||||
SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples)
|
||||
{
|
||||
int numRows = pvdgTable.numRows();
|
||||
|
||||
// Copy data
|
||||
const std::vector<double> &press = pvdgTable.getPressureColumn();
|
||||
const std::vector<double> &B = pvdgTable.getFormationFactorColumn();
|
||||
const std::vector<double> &visc = pvdgTable.getViscosityColumn();
|
||||
|
||||
std::vector<double> B_inv(numRows);
|
||||
for (int i = 0; i < numRows; ++i) {
|
||||
B_inv[i] = 1.0 / B[i];
|
||||
}
|
||||
|
||||
buildUniformMonotoneTable(press, B_inv, samples, b_);
|
||||
buildUniformMonotoneTable(press, visc, samples, viscosity_);
|
||||
}
|
||||
|
||||
// Destructor
|
||||
SinglePvtDeadSpline::~SinglePvtDeadSpline()
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
|
||||
void SinglePvtDeadSpline::mu(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_mu) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_mu[i] = viscosity_(p[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtDeadSpline::mu(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_mu[i] = viscosity_(p[i]);
|
||||
output_dmudp[i] = viscosity_.derivative(p[i]);
|
||||
}
|
||||
std::fill(output_dmudr, output_dmudr + n, 0.0);
|
||||
|
||||
}
|
||||
|
||||
void SinglePvtDeadSpline::mu(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
const PhasePresence* /*cond*/,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_mu[i] = viscosity_(p[i]);
|
||||
output_dmudp[i] = viscosity_.derivative(p[i]);
|
||||
}
|
||||
std::fill(output_dmudr, output_dmudr + n, 0.0);
|
||||
|
||||
}
|
||||
|
||||
void SinglePvtDeadSpline::B(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_B) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_B[i] = 1.0/b_(p[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtDeadSpline::dBdp(const int n,
|
||||
const double* p,
|
||||
const double* /*z*/,
|
||||
double* output_B,
|
||||
double* output_dBdp) const
|
||||
{
|
||||
B(n, p, 0, output_B);
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
double Bg = output_B[i];
|
||||
output_dBdp[i] = -Bg*Bg*b_.derivative(p[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtDeadSpline::b(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_b[i] = b_(p[i]);
|
||||
output_dbdp[i] = b_.derivative(p[i]);
|
||||
|
||||
}
|
||||
std::fill(output_dbdr, output_dbdr + n, 0.0);
|
||||
|
||||
}
|
||||
|
||||
void SinglePvtDeadSpline::b(const int n,
|
||||
const double* p,
|
||||
const double* /*r*/,
|
||||
const PhasePresence* /*cond*/,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_b[i] = b_(p[i]);
|
||||
output_dbdp[i] = b_.derivative(p[i]);
|
||||
|
||||
}
|
||||
std::fill(output_dbdr, output_dbdr + n, 0.0);
|
||||
|
||||
}
|
||||
|
||||
void SinglePvtDeadSpline::rsSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rsSat,
|
||||
double* output_drsSatdp) const
|
||||
{
|
||||
std::fill(output_rsSat, output_rsSat + n, 0.0);
|
||||
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
void SinglePvtDeadSpline::rvSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rvSat,
|
||||
double* output_drvSatdp) const
|
||||
{
|
||||
std::fill(output_rvSat, output_rvSat + n, 0.0);
|
||||
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
void SinglePvtDeadSpline::R(const int n,
|
||||
const double* /*p*/,
|
||||
const double* /*z*/,
|
||||
double* output_R) const
|
||||
{
|
||||
std::fill(output_R, output_R + n, 0.0);
|
||||
}
|
||||
|
||||
void SinglePvtDeadSpline::dRdp(const int n,
|
||||
const double* /*p*/,
|
||||
const double* /*z*/,
|
||||
double* output_R,
|
||||
double* output_dRdp) const
|
||||
{
|
||||
std::fill(output_R, output_R + n, 0.0);
|
||||
std::fill(output_dRdp, output_dRdp + n, 0.0);
|
||||
}
|
||||
|
||||
}
|
||||
@ -1,48 +0,0 @@
|
||||
/*
|
||||
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
|
||||
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
SinglePvtInterface::SinglePvtInterface()
|
||||
: num_phases_(MaxNumPhases)
|
||||
{
|
||||
for (int i = 0; i < MaxNumPhases; ++i) {
|
||||
phase_pos_[i] = i;
|
||||
}
|
||||
}
|
||||
|
||||
SinglePvtInterface::~SinglePvtInterface()
|
||||
{
|
||||
}
|
||||
|
||||
void SinglePvtInterface::setPhaseConfiguration(const int num_phases, const int* phase_pos)
|
||||
{
|
||||
num_phases_ = num_phases;
|
||||
for (int i = 0; i < MaxNumPhases; ++i) {
|
||||
phase_pos_[i] = phase_pos[i];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
@ -1,492 +0,0 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: MiscibilityLiveGas.cpp
|
||||
//
|
||||
// Created: Wed Feb 10 09:21:53 2010
|
||||
//
|
||||
// Author: Bjørn Spjelkavik <bsp@sintef.no>
|
||||
//
|
||||
// Revision: $Id$
|
||||
//
|
||||
//===========================================================================
|
||||
/*
|
||||
Copyright 2010 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/pvt/SinglePvtLiveGas.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/linearInterpolation.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
using Opm::linearInterpolation;
|
||||
using Opm::linearInterpolationDerivative;
|
||||
|
||||
|
||||
//------------------------------------------------------------------------
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
SinglePvtLiveGas::SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable)
|
||||
{
|
||||
// GAS, PVTG
|
||||
saturated_gas_table_.resize(4);
|
||||
saturated_gas_table_[0] = pvtgTable.getOuterTable()->getPressureColumn();
|
||||
saturated_gas_table_[1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn();
|
||||
saturated_gas_table_[2] = pvtgTable.getOuterTable()->getGasViscosityColumn();
|
||||
saturated_gas_table_[3] = pvtgTable.getOuterTable()->getOilSolubilityColumn();
|
||||
|
||||
int sz = pvtgTable.getOuterTable()->numRows();
|
||||
undersat_gas_tables_.resize(sz);
|
||||
for (int i=0; i<sz; ++i) {
|
||||
const auto &undersatTable = *pvtgTable.getInnerTable(i);
|
||||
|
||||
undersat_gas_tables_[i].resize(3);
|
||||
undersat_gas_tables_[i][0] = undersatTable.getOilSolubilityColumn();
|
||||
undersat_gas_tables_[i][1] = undersatTable.getGasFormationFactorColumn();
|
||||
undersat_gas_tables_[i][2] = undersatTable.getGasViscosityColumn();
|
||||
}
|
||||
|
||||
// Bg -> 1/Bg
|
||||
for (int i=0; i<sz; ++i) {
|
||||
saturated_gas_table_[1][i] = 1.0/saturated_gas_table_[1][i];
|
||||
for (size_t j=0; j<undersat_gas_tables_[i][1].size(); ++j) {
|
||||
undersat_gas_tables_[i][1][j] = 1.0/undersat_gas_tables_[i][1][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Destructor
|
||||
SinglePvtLiveGas::~SinglePvtLiveGas()
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
void SinglePvtLiveGas::mu(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_mu) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_mu[i] = miscible_gas(p[i], z + num_phases_*i, 2, false);
|
||||
}
|
||||
}
|
||||
|
||||
/// Viscosity and its derivatives as a function of p and r.
|
||||
void SinglePvtLiveGas::mu(const int /*n*/,
|
||||
const double* /*p*/,
|
||||
const double* /*r*/,
|
||||
double* /*output_mu*/,
|
||||
double* /*output_dmudp*/,
|
||||
double* /*output_dmudr*/) const
|
||||
{
|
||||
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
|
||||
}
|
||||
|
||||
/// Viscosity and its derivatives as a function of p and r.
|
||||
void SinglePvtLiveGas::mu(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
const PhasePresence* cond,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const PhasePresence& cnd = cond[i];
|
||||
output_mu[i] = miscible_gas(p[i], r[i], cnd,2, 0);
|
||||
output_dmudp[i] = miscible_gas(p[i], r[i], cnd, 2, 1);
|
||||
output_dmudr[i] = miscible_gas(p[i], r[i], cnd, 2, 2);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/// Formation volume factor as a function of p and z.
|
||||
void SinglePvtLiveGas::B(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_B) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_B[i] = evalB(p[i], z + num_phases_*i);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/// Formation volume factor and p-derivative as functions of p and z.
|
||||
void SinglePvtLiveGas::dBdp(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_B,
|
||||
double* output_dBdp) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
evalBDeriv(p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
|
||||
}
|
||||
}
|
||||
|
||||
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
|
||||
void SinglePvtLiveGas::b(const int /*n*/,
|
||||
const double* /*p*/,
|
||||
const double* /*r*/,
|
||||
double* /*output_b*/,
|
||||
double* /*output_dbdp*/,
|
||||
double* /*output_dbdr*/) const
|
||||
|
||||
{
|
||||
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
|
||||
}
|
||||
|
||||
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
|
||||
void SinglePvtLiveGas::b(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
const PhasePresence* cond,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const PhasePresence& cnd = cond[i];
|
||||
|
||||
output_b[i] = miscible_gas(p[i], r[i], cnd, 1, 0);
|
||||
output_dbdp[i] = miscible_gas(p[i], r[i], cnd, 1, 1);
|
||||
output_dbdr[i] = miscible_gas(p[i], r[i], cnd, 1, 2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
/// Gas resolution and its derivatives at bublepoint as a function of p.
|
||||
void SinglePvtLiveGas::rvSat(const int n,
|
||||
const double* p,
|
||||
double* output_rvSat,
|
||||
double* output_drvSatdp) const
|
||||
{
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_rvSat[i] = linearInterpolation(saturated_gas_table_[0],
|
||||
saturated_gas_table_[3],p[i]);
|
||||
output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[0],
|
||||
saturated_gas_table_[3],p[i]);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtLiveGas::rsSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rsSat,
|
||||
double* output_drsSatdp) const
|
||||
{
|
||||
std::fill(output_rsSat, output_rsSat + n, 0.0);
|
||||
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
/// Solution factor as a function of p and z.
|
||||
void SinglePvtLiveGas::R(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_R) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_R[i] = evalR(p[i], z + num_phases_*i);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/// Solution factor and p-derivative as functions of p and z.
|
||||
void SinglePvtLiveGas::dRdp(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_R,
|
||||
double* output_dRdp) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
evalRDeriv(p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ---- Private methods ----
|
||||
|
||||
double SinglePvtLiveGas::evalB(const double press, const double* surfvol) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Vapour]] == 0.0) {
|
||||
// To handle no-gas case.
|
||||
return 1.0;
|
||||
}
|
||||
return 1.0/miscible_gas(press, surfvol, 1, false);
|
||||
}
|
||||
|
||||
void SinglePvtLiveGas::evalBDeriv(const double press, const double* surfvol,
|
||||
double& Bval, double& dBdpval) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Vapour]] == 0.0) {
|
||||
// To handle no-gas case.
|
||||
Bval = 1.0;
|
||||
dBdpval = 0.0;
|
||||
return;
|
||||
}
|
||||
Bval = evalB(press, surfvol);
|
||||
dBdpval = -Bval*Bval*miscible_gas(press, surfvol, 1, true);
|
||||
}
|
||||
|
||||
double SinglePvtLiveGas::evalR(const double press, const double* surfvol) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Liquid]] == 0.0) {
|
||||
// To handle no-gas case.
|
||||
return 0.0;
|
||||
}
|
||||
double satR = linearInterpolation(saturated_gas_table_[0],
|
||||
saturated_gas_table_[3], press);
|
||||
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
|
||||
if (satR < maxR ) {
|
||||
// Saturated case
|
||||
return satR;
|
||||
} else {
|
||||
// Undersaturated case
|
||||
return maxR;
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtLiveGas::evalRDeriv(const double press, const double* surfvol,
|
||||
double& Rval, double& dRdpval) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Liquid]] == 0.0) {
|
||||
// To handle no-gas case.
|
||||
Rval = 0.0;
|
||||
dRdpval = 0.0;
|
||||
return;
|
||||
}
|
||||
double satR = linearInterpolation(saturated_gas_table_[0],
|
||||
saturated_gas_table_[3], press);
|
||||
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
|
||||
if (satR < maxR ) {
|
||||
// Saturated case
|
||||
Rval = satR;
|
||||
dRdpval = linearInterpolationDerivative(saturated_gas_table_[0],
|
||||
saturated_gas_table_[3],
|
||||
press);
|
||||
} else {
|
||||
// Undersaturated case
|
||||
Rval = maxR;
|
||||
dRdpval = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
double SinglePvtLiveGas::miscible_gas(const double press,
|
||||
const double* surfvol,
|
||||
const int item,
|
||||
const bool deriv) const
|
||||
{
|
||||
int section;
|
||||
double Rval = linearInterpolation(saturated_gas_table_[0],
|
||||
saturated_gas_table_[3], press,
|
||||
section);
|
||||
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
|
||||
if (deriv) {
|
||||
if (Rval < maxR ) { // Saturated case
|
||||
return linearInterpolationDerivative(saturated_gas_table_[0],
|
||||
saturated_gas_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = section;
|
||||
if (undersat_gas_tables_[is][0].size() < 2) {
|
||||
double val = (saturated_gas_table_[item][is+1]
|
||||
- saturated_gas_table_[item][is]) /
|
||||
(saturated_gas_table_[0][is+1] -
|
||||
saturated_gas_table_[0][is]);
|
||||
return val;
|
||||
}
|
||||
double val1 =
|
||||
linearInterpolation(undersat_gas_tables_[is][0],
|
||||
undersat_gas_tables_[is][item],
|
||||
maxR);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_gas_tables_[is+1][0],
|
||||
undersat_gas_tables_[is+1][item],
|
||||
maxR);
|
||||
double val = (val2 - val1)/
|
||||
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
|
||||
return val;
|
||||
}
|
||||
} else {
|
||||
if (Rval < maxR ) { // Saturated case
|
||||
return linearInterpolation(saturated_gas_table_[0],
|
||||
saturated_gas_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = section;
|
||||
// Extrapolate from first table section
|
||||
if (is == 0 && press < saturated_gas_table_[0][0]) {
|
||||
return linearInterpolation(undersat_gas_tables_[0][0],
|
||||
undersat_gas_tables_[0][item],
|
||||
maxR);
|
||||
}
|
||||
|
||||
// Extrapolate from last table section
|
||||
int ltp = saturated_gas_table_[0].size() - 1;
|
||||
if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
|
||||
return linearInterpolation(undersat_gas_tables_[ltp][0],
|
||||
undersat_gas_tables_[ltp][item],
|
||||
maxR);
|
||||
}
|
||||
|
||||
// Interpolate between table sections
|
||||
double w = (press - saturated_gas_table_[0][is]) /
|
||||
(saturated_gas_table_[0][is+1] -
|
||||
saturated_gas_table_[0][is]);
|
||||
if (undersat_gas_tables_[is][0].size() < 2) {
|
||||
double val = saturated_gas_table_[item][is] +
|
||||
w*(saturated_gas_table_[item][is+1] -
|
||||
saturated_gas_table_[item][is]);
|
||||
return val;
|
||||
}
|
||||
double val1 =
|
||||
linearInterpolation(undersat_gas_tables_[is][0],
|
||||
undersat_gas_tables_[is][item],
|
||||
maxR);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_gas_tables_[is+1][0],
|
||||
undersat_gas_tables_[is+1][item],
|
||||
maxR);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double SinglePvtLiveGas::miscible_gas(const double press,
|
||||
const double r,
|
||||
const PhasePresence& cond,
|
||||
const int item,
|
||||
const int deriv) const
|
||||
{
|
||||
const bool isSat = cond.hasFreeOil();
|
||||
|
||||
// Derivative w.r.t p
|
||||
if (deriv == 1) {
|
||||
if (isSat) { // Saturated case
|
||||
return linearInterpolationDerivative(saturated_gas_table_[0],
|
||||
saturated_gas_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_gas_table_[0], press);
|
||||
if (undersat_gas_tables_[is][0].size() < 2) {
|
||||
double val = (saturated_gas_table_[item][is+1]
|
||||
- saturated_gas_table_[item][is]) /
|
||||
(saturated_gas_table_[0][is+1] -
|
||||
saturated_gas_table_[0][is]);
|
||||
return val;
|
||||
}
|
||||
double val1 =
|
||||
linearInterpolation(undersat_gas_tables_[is][0],
|
||||
undersat_gas_tables_[is][item],
|
||||
r);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_gas_tables_[is+1][0],
|
||||
undersat_gas_tables_[is+1][item],
|
||||
r);
|
||||
double val = (val2 - val1)/
|
||||
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
|
||||
return val;
|
||||
}
|
||||
} else if (deriv == 2){
|
||||
if (isSat) {
|
||||
return 0;
|
||||
} else {
|
||||
int is = tableIndex(saturated_gas_table_[0], press);
|
||||
double w = (press - saturated_gas_table_[0][is]) /
|
||||
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
|
||||
assert(undersat_gas_tables_[is][0].size() >= 2);
|
||||
assert(undersat_gas_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolationDerivative(undersat_gas_tables_[is][0],
|
||||
undersat_gas_tables_[is][item],
|
||||
r);
|
||||
double val2 =
|
||||
linearInterpolationDerivative(undersat_gas_tables_[is+1][0],
|
||||
undersat_gas_tables_[is+1][item],
|
||||
r);
|
||||
|
||||
double val = val1 + w * (val2 - val1);
|
||||
return val;
|
||||
|
||||
}
|
||||
} else {
|
||||
if (isSat) { // Saturated case
|
||||
return linearInterpolation(saturated_gas_table_[0],
|
||||
saturated_gas_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_gas_table_[0], press);
|
||||
// Extrapolate from first table section
|
||||
if (is == 0 && press < saturated_gas_table_[0][0]) {
|
||||
return linearInterpolation(undersat_gas_tables_[0][0],
|
||||
undersat_gas_tables_[0][item],
|
||||
r);
|
||||
}
|
||||
|
||||
// Extrapolate from last table section
|
||||
//int ltp = saturated_gas_table_[0].size() - 1;
|
||||
//if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
|
||||
// return linearInterpolation(undersat_gas_tables_[ltp][0],
|
||||
// undersat_gas_tables_[ltp][item],
|
||||
// r);
|
||||
//}
|
||||
|
||||
// Interpolate between table sections
|
||||
double w = (press - saturated_gas_table_[0][is]) /
|
||||
(saturated_gas_table_[0][is+1] -
|
||||
saturated_gas_table_[0][is]);
|
||||
if (undersat_gas_tables_[is][0].size() < 2) {
|
||||
double val = saturated_gas_table_[item][is] +
|
||||
w*(saturated_gas_table_[item][is+1] -
|
||||
saturated_gas_table_[item][is]);
|
||||
return val;
|
||||
}
|
||||
double val1 =
|
||||
linearInterpolation(undersat_gas_tables_[is][0],
|
||||
undersat_gas_tables_[is][item],
|
||||
r);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_gas_tables_[is+1][0],
|
||||
undersat_gas_tables_[is+1][item],
|
||||
r);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
@ -1,551 +0,0 @@
|
||||
/*
|
||||
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/linearInterpolation.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
using Opm::linearInterpolation;
|
||||
using Opm::linearInterpolationDerivative;
|
||||
using Opm::tableIndex;
|
||||
|
||||
//------------------------------------------------------------------------
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
|
||||
SinglePvtLiveOil::SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable)
|
||||
{
|
||||
const auto saturatedPvto = pvtoTable.getOuterTable();
|
||||
|
||||
// OIL, PVTO
|
||||
saturated_oil_table_.resize(4);
|
||||
const int sz = saturatedPvto->numRows();
|
||||
for (int k=0; k<4; ++k) {
|
||||
saturated_oil_table_[k].resize(sz);
|
||||
}
|
||||
for (int i=0; i<sz; ++i) {
|
||||
saturated_oil_table_[0][i] = saturatedPvto->getPressureColumn()[i]; // p
|
||||
saturated_oil_table_[1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo
|
||||
saturated_oil_table_[2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o
|
||||
saturated_oil_table_[3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs
|
||||
}
|
||||
|
||||
undersat_oil_tables_.resize(sz);
|
||||
for (int i=0; i<sz; ++i) {
|
||||
const auto undersaturatedPvto = pvtoTable.getInnerTable(i);
|
||||
|
||||
undersat_oil_tables_[i].resize(3);
|
||||
int tsize = undersaturatedPvto->numRows();
|
||||
undersat_oil_tables_[i][0].resize(tsize);
|
||||
undersat_oil_tables_[i][1].resize(tsize);
|
||||
undersat_oil_tables_[i][2].resize(tsize);
|
||||
for (int j=0; j<tsize; ++j) {
|
||||
undersat_oil_tables_[i][0][j] = undersaturatedPvto->getPressureColumn()[j]; // p
|
||||
undersat_oil_tables_[i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j]; // 1/Bo
|
||||
undersat_oil_tables_[i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j]; // mu_o
|
||||
}
|
||||
}
|
||||
|
||||
// Complete undersaturated tables by extrapolating from existing data
|
||||
// as is done in Eclipse and Mrst
|
||||
int iNext = -1;
|
||||
for (int i=0; i<sz; ++i) {
|
||||
// Skip records already containing undersaturated data
|
||||
if (undersat_oil_tables_[i][0].size() > 1) {
|
||||
continue;
|
||||
}
|
||||
// Look ahead for next record containing undersaturated data
|
||||
if (iNext < i) {
|
||||
iNext = i+1;
|
||||
while (iNext<sz && undersat_oil_tables_[iNext][0].size() < 2) {
|
||||
++iNext;
|
||||
}
|
||||
if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
|
||||
}
|
||||
// Add undersaturated data to current record while maintaining compressibility and viscosibility
|
||||
typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
|
||||
for (sz_t j=1; j<undersat_oil_tables_[iNext][0].size(); ++j) {
|
||||
double diffPressure = undersat_oil_tables_[iNext][0][j]-undersat_oil_tables_[iNext][0][j-1];
|
||||
double pressure = undersat_oil_tables_[i][0].back()+diffPressure;
|
||||
undersat_oil_tables_[i][0].push_back(pressure);
|
||||
double compr = (1.0/undersat_oil_tables_[iNext][1][j]-1.0/undersat_oil_tables_[iNext][1][j-1])
|
||||
/ (0.5*(1.0/undersat_oil_tables_[iNext][1][j]+1.0/undersat_oil_tables_[iNext][1][j-1]));
|
||||
double B = (1.0/undersat_oil_tables_[i][1].back())*(1.0+0.5*compr)/(1.0-0.5*compr);
|
||||
undersat_oil_tables_[i][1].push_back(1.0/B);
|
||||
double visc = (undersat_oil_tables_[iNext][2][j]-undersat_oil_tables_[iNext][2][j-1])
|
||||
/ (0.5*(undersat_oil_tables_[iNext][2][j]+undersat_oil_tables_[iNext][2][j-1]));
|
||||
double mu = (undersat_oil_tables_[i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc);
|
||||
undersat_oil_tables_[i][2].push_back(mu);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Destructor.
|
||||
SinglePvtLiveOil::~SinglePvtLiveOil()
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
/// Viscosity as a function of p and z.
|
||||
void SinglePvtLiveOil::mu(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_mu) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_mu[i] = miscible_oil(p[i], z + num_phases_*i, 2, false);
|
||||
}
|
||||
}
|
||||
|
||||
/// Viscosity and its derivatives as a function of p and r.
|
||||
void SinglePvtLiveOil::mu(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_mu[i] = miscible_oil(p[i], r[i], 2, 0);
|
||||
output_dmudp[i] = miscible_oil(p[i], r[i], 2, 1);
|
||||
output_dmudr[i] = miscible_oil(p[i], r[i], 2, 2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
/// Viscosity and its derivatives as a function of p and r.
|
||||
void SinglePvtLiveOil::mu(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
const PhasePresence* cond,
|
||||
double* output_mu,
|
||||
double* output_dmudp,
|
||||
double* output_dmudr) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const PhasePresence& cnd = cond[i];
|
||||
|
||||
output_mu[i] = miscible_oil(p[i], r[i], cnd, 2, 0);
|
||||
output_dmudp[i] = miscible_oil(p[i], r[i], cnd, 2, 1);
|
||||
output_dmudr[i] = miscible_oil(p[i], r[i], cnd, 2, 2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// Formation volume factor as a function of p and z.
|
||||
void SinglePvtLiveOil::B(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_B) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_B[i] = evalB(p[i], z + num_phases_*i);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/// Formation volume factor and p-derivative as functions of p and z.
|
||||
void SinglePvtLiveOil::dBdp(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_B,
|
||||
double* output_dBdp) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
evalBDeriv(p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtLiveOil::b(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_b[i] = miscible_oil(p[i], r[i], 1, 0);
|
||||
output_dbdp[i] = miscible_oil(p[i], r[i], 1, 1);
|
||||
output_dbdr[i] = miscible_oil(p[i], r[i], 1, 2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtLiveOil::b(const int n,
|
||||
const double* p,
|
||||
const double* r,
|
||||
const PhasePresence* cond,
|
||||
double* output_b,
|
||||
double* output_dbdp,
|
||||
double* output_dbdr) const
|
||||
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const PhasePresence& cnd = cond[i];
|
||||
|
||||
output_b[i] = miscible_oil(p[i], r[i], cnd, 1, 0);
|
||||
output_dbdp[i] = miscible_oil(p[i], r[i], cnd, 1, 1);
|
||||
output_dbdr[i] = miscible_oil(p[i], r[i], cnd, 1, 2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtLiveOil::rsSat(const int n,
|
||||
const double* p,
|
||||
double* output_rsSat,
|
||||
double* output_drsSatdp) const
|
||||
{
|
||||
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_rsSat[i] = linearInterpolation(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3],p[i]);
|
||||
output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3],p[i]);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtLiveOil::rvSat(const int n,
|
||||
const double* /*p*/,
|
||||
double* output_rvSat,
|
||||
double* output_drvSatdp) const
|
||||
{
|
||||
std::fill(output_rvSat, output_rvSat + n, 0.0);
|
||||
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
|
||||
}
|
||||
|
||||
/// Solution factor as a function of p and z.
|
||||
void SinglePvtLiveOil::R(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_R) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
output_R[i] = evalR(p[i], z + num_phases_*i);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/// Solution factor and p-derivative as functions of p and z.
|
||||
void SinglePvtLiveOil::dRdp(const int n,
|
||||
const double* p,
|
||||
const double* z,
|
||||
double* output_R,
|
||||
double* output_dRdp) const
|
||||
{
|
||||
// #pragma omp parallel for
|
||||
for (int i = 0; i < n; ++i) {
|
||||
evalRDeriv(p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
// ---- Private methods ----
|
||||
|
||||
double SinglePvtLiveOil::evalB(double press, const double* surfvol) const
|
||||
{
|
||||
// if (surfvol[phase_pos_[Liquid]] == 0.0) return 1.0; // To handle no-oil case.
|
||||
return 1.0/miscible_oil(press, surfvol, 1, false);
|
||||
}
|
||||
|
||||
|
||||
void SinglePvtLiveOil::evalBDeriv(const double press, const double* surfvol,
|
||||
double& Bval, double& dBdpval) const
|
||||
{
|
||||
Bval = evalB(press, surfvol);
|
||||
dBdpval = -Bval*Bval*miscible_oil(press, surfvol, 1, true);
|
||||
}
|
||||
|
||||
double SinglePvtLiveOil::evalR(double press, const double* surfvol) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Vapour]] == 0.0) {
|
||||
return 0.0;
|
||||
}
|
||||
double Rval = linearInterpolation(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3], press);
|
||||
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
|
||||
if (Rval < maxR ) { // Saturated case
|
||||
return Rval;
|
||||
} else {
|
||||
return maxR; // Undersaturated case
|
||||
}
|
||||
}
|
||||
|
||||
void SinglePvtLiveOil::evalRDeriv(const double press, const double* surfvol,
|
||||
double& Rval, double& dRdpval) const
|
||||
{
|
||||
if (surfvol[phase_pos_[Vapour]] == 0.0) {
|
||||
Rval = 0.0;
|
||||
dRdpval = 0.0;
|
||||
return;
|
||||
}
|
||||
Rval = linearInterpolation(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3], press);
|
||||
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
|
||||
if (Rval < maxR ) {
|
||||
// Saturated case
|
||||
dRdpval = linearInterpolationDerivative(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3],
|
||||
press);
|
||||
} else {
|
||||
// Undersaturated case
|
||||
Rval = maxR;
|
||||
dRdpval = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double SinglePvtLiveOil::miscible_oil(const double press,
|
||||
const double* surfvol,
|
||||
const int item,
|
||||
const bool deriv) const
|
||||
{
|
||||
int section;
|
||||
double Rval = linearInterpolation(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3],
|
||||
press, section);
|
||||
double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
|
||||
if (deriv) {
|
||||
if (Rval < maxR ) { // Saturated case
|
||||
return linearInterpolationDerivative(saturated_oil_table_[0],
|
||||
saturated_oil_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[3], maxR);
|
||||
double w = (maxR - saturated_oil_table_[3][is]) /
|
||||
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
|
||||
assert(undersat_oil_tables_[is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[is][0],
|
||||
undersat_oil_tables_[is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
|
||||
undersat_oil_tables_[is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
} else {
|
||||
if (Rval < maxR ) { // Saturated case
|
||||
return linearInterpolation(saturated_oil_table_[0],
|
||||
saturated_oil_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
// Interpolate between table sections
|
||||
int is = tableIndex(saturated_oil_table_[3], maxR);
|
||||
double w = (maxR - saturated_oil_table_[3][is]) /
|
||||
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
|
||||
assert(undersat_oil_tables_[is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolation(undersat_oil_tables_[is][0],
|
||||
undersat_oil_tables_[is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_oil_tables_[is+1][0],
|
||||
undersat_oil_tables_[is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double SinglePvtLiveOil::miscible_oil(const double press,
|
||||
const double r,
|
||||
const int item,
|
||||
const int deriv) const
|
||||
{
|
||||
int section;
|
||||
double Rval = linearInterpolation(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3],
|
||||
press, section);
|
||||
// derivative with respect to frist component (pressure)
|
||||
if (deriv == 1) {
|
||||
if (Rval <= r ) { // Saturated case
|
||||
return linearInterpolationDerivative(saturated_oil_table_[0],
|
||||
saturated_oil_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[3], r);
|
||||
double w = (r - saturated_oil_table_[3][is]) /
|
||||
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
|
||||
assert(undersat_oil_tables_[is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[is][0],
|
||||
undersat_oil_tables_[is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
|
||||
undersat_oil_tables_[is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
// derivative with respect to second component (r)
|
||||
} else if (deriv == 2) {
|
||||
if (Rval <= r ) { // Saturated case
|
||||
return 0;
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[3], r);
|
||||
assert(undersat_oil_tables_[is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolation(undersat_oil_tables_[is][0],
|
||||
undersat_oil_tables_[is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_oil_tables_[is+1][0],
|
||||
undersat_oil_tables_[is+1][item],
|
||||
press);
|
||||
|
||||
double val = (val2 - val1)/(saturated_oil_table_[3][is+1]-saturated_oil_table_[3][is]);
|
||||
return val;
|
||||
}
|
||||
|
||||
|
||||
} else {
|
||||
if (Rval <= r ) { // Saturated case
|
||||
return linearInterpolation(saturated_oil_table_[0],
|
||||
saturated_oil_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
// Interpolate between table sections
|
||||
int is = tableIndex(saturated_oil_table_[3], r);
|
||||
double w = (r - saturated_oil_table_[3][is]) /
|
||||
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
|
||||
assert(undersat_oil_tables_[is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolation(undersat_oil_tables_[is][0],
|
||||
undersat_oil_tables_[is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_oil_tables_[is+1][0],
|
||||
undersat_oil_tables_[is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double SinglePvtLiveOil::miscible_oil(const double press,
|
||||
const double r,
|
||||
const PhasePresence& cond,
|
||||
const int item,
|
||||
const int deriv) const
|
||||
{
|
||||
const bool isSat = cond.hasFreeGas();
|
||||
|
||||
// derivative with respect to frist component (pressure)
|
||||
if (deriv == 1) {
|
||||
if (isSat) { // Saturated case
|
||||
return linearInterpolationDerivative(saturated_oil_table_[0],
|
||||
saturated_oil_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[3], r);
|
||||
double w = (r - saturated_oil_table_[3][is]) /
|
||||
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
|
||||
assert(undersat_oil_tables_[is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[is][0],
|
||||
undersat_oil_tables_[is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
|
||||
undersat_oil_tables_[is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
// derivative with respect to second component (r)
|
||||
} else if (deriv == 2) {
|
||||
if (isSat) { // Saturated case
|
||||
return 0;
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[3], r);
|
||||
assert(undersat_oil_tables_[is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolation(undersat_oil_tables_[is][0],
|
||||
undersat_oil_tables_[is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_oil_tables_[is+1][0],
|
||||
undersat_oil_tables_[is+1][item],
|
||||
press);
|
||||
|
||||
double val = (val2 - val1)/(saturated_oil_table_[3][is+1]-saturated_oil_table_[3][is]);
|
||||
return val;
|
||||
}
|
||||
|
||||
|
||||
} else {
|
||||
if (isSat) { // Saturated case
|
||||
return linearInterpolation(saturated_oil_table_[0],
|
||||
saturated_oil_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
// Interpolate between table sections
|
||||
int is = tableIndex(saturated_oil_table_[3], r);
|
||||
double w = (r - saturated_oil_table_[3][is]) /
|
||||
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
|
||||
assert(undersat_oil_tables_[is][0].size() >= 2);
|
||||
assert(undersat_oil_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolation(undersat_oil_tables_[is][0],
|
||||
undersat_oil_tables_[is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolation(undersat_oil_tables_[is+1][0],
|
||||
undersat_oil_tables_[is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
@ -279,7 +279,7 @@ namespace Opm
|
||||
if (rec[i].live_oil_table_index > 0) {
|
||||
if (deck->hasKeyword("RSVD") &&
|
||||
size_t(rec[i].live_oil_table_index) <= deck->getKeyword("RSVD")->size()) {
|
||||
Opm::SingleRecordTable rsvd(deck->getKeyword("RSVD"),std::vector<std::string>{"vd", "rs"},rec[i].live_oil_table_index-1);
|
||||
Opm::SingleRecordTable rsvd(deck->getKeyword("RSVD"),std::vector<std::string>{"vd", "rs"},rec[i].live_oil_table_index-1);
|
||||
std::vector<double> vd(rsvd.getColumn("vd"));
|
||||
std::vector<double> rs(rsvd.getColumn("rs"));
|
||||
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, vd, rs));
|
||||
|
||||
@ -9,6 +9,19 @@ WATER
|
||||
|
||||
FIELD
|
||||
|
||||
-- tests for the PVT functions need a grid because the OPM-API for the
|
||||
-- PVT functions assumes the presence of compressed grid cells,
|
||||
-- i.e. we need to be able to map from compressed to logical cartesian
|
||||
-- cell indices to find out the PVT region of a cell
|
||||
DIMENS
|
||||
3 3 3 /
|
||||
DXV
|
||||
1 2 3 /
|
||||
DYV
|
||||
1 2 3 /
|
||||
DZV
|
||||
1 2 3 /
|
||||
|
||||
PVTO
|
||||
-- Rs Pbub Bo Vo
|
||||
.0 14.7 1.0000 1.20 /
|
||||
|
||||
@ -1,15 +1,17 @@
|
||||
#include <config.h>
|
||||
|
||||
#include <opm/core/grid/GridManager.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtConstCompr.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtDead.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
|
||||
#include <opm/core/props/pvt/PvtConstCompr.hpp>
|
||||
#include <opm/core/props/pvt/PvtDead.hpp>
|
||||
#include <opm/core/props/pvt/PvtDeadSpline.hpp>
|
||||
#include <opm/core/props/pvt/PvtLiveOil.hpp>
|
||||
#include <opm/core/props/pvt/PvtLiveGas.hpp>
|
||||
#include <opm/core/props/phaseUsageFromDeck.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
|
||||
#include <opm/core/utility/Units.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/grid.h>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
|
||||
#include <opm/parser/eclipse/Utility/PvtgTable.hpp>
|
||||
@ -36,63 +38,78 @@ using namespace Opm;
|
||||
using namespace std;
|
||||
|
||||
|
||||
std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr deck,PhaseUsage phase_usage_){
|
||||
|
||||
std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, PhaseUsage phase_usage_){
|
||||
Opm::GridManager grid(deck);
|
||||
const UnstructuredGrid* cgrid = grid.c_grid();
|
||||
|
||||
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
|
||||
int samples = 0;
|
||||
|
||||
// first, calculate the PVT table index for each compressed
|
||||
// cell. This array is required to construct the PVT classes
|
||||
// below.
|
||||
std::vector<int> pvtTableIdx;
|
||||
Opm::extractPvtTableIndex(pvtTableIdx,
|
||||
deck,
|
||||
cgrid->number_of_cells,
|
||||
cgrid->global_cell);
|
||||
|
||||
std::vector<std::shared_ptr<SinglePvtInterface> > props_;
|
||||
std::vector<std::shared_ptr<PvtInterface> > props_;
|
||||
// Set the properties.
|
||||
props_.resize(phase_usage_.num_phases);
|
||||
|
||||
// Water PVT
|
||||
if (phase_usage_.phase_used[Aqua]) {
|
||||
if (deck->hasKeyword("PVTW")) {
|
||||
Opm::PvtwTable pvtwTable(deck->getKeyword("PVTW"));
|
||||
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable));
|
||||
std::shared_ptr<PvtConstCompr> pvtw(new PvtConstCompr);
|
||||
pvtw->initFromWater(deck->getKeyword("PVTW"), pvtTableIdx);
|
||||
|
||||
props_[phase_usage_.phase_pos[Aqua]] = pvtw;
|
||||
} else {
|
||||
// Eclipse 100 default.
|
||||
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
|
||||
props_[phase_usage_.phase_pos[Aqua]].reset(new PvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
|
||||
}
|
||||
}
|
||||
|
||||
// Oil PVT
|
||||
if (phase_usage_.phase_used[Liquid]) {
|
||||
if (deck->hasKeyword("PVDO")) {
|
||||
Opm::PvdoTable pvdoTable(deck->getKeyword("PVDO"));
|
||||
Opm::DeckKeywordConstPtr pvdoKeyword(deck->getKeyword("PVDO"));
|
||||
if (samples > 0) {
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
|
||||
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
|
||||
splinePvt->initFromOil(pvdoKeyword, pvtTableIdx, samples);
|
||||
props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
|
||||
} else {
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable));
|
||||
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
|
||||
deadPvt->initFromOil(pvdoKeyword, pvtTableIdx);
|
||||
props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
|
||||
}
|
||||
} else if (deck->hasKeyword("PVTO")) {
|
||||
Opm::PvtoTable pvtoTable(deck->getKeyword("PVTO"), /*tableIdx=*/0);
|
||||
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable));
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"), pvtTableIdx));
|
||||
} else if (deck->hasKeyword("PVCDO")) {
|
||||
Opm::PvcdoTable pvcdoTable(deck->getKeyword("PVCDO"), /*tableIdx=*/0);
|
||||
std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
|
||||
pvcdo->initFromOil(deck->getKeyword("PVCDO"), pvtTableIdx);
|
||||
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable));
|
||||
props_[phase_usage_.phase_pos[Liquid]] = pvcdo;
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
|
||||
OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n");
|
||||
}
|
||||
}
|
||||
// Gas PVT
|
||||
if (phase_usage_.phase_used[Vapour]) {
|
||||
if (deck->hasKeyword("PVDG")) {
|
||||
Opm::PvdgTable pvdgTable(deck->getKeyword("PVDG"));
|
||||
|
||||
Opm::DeckKeywordConstPtr pvdgKeyword(deck->getKeyword("PVDG"));
|
||||
if (samples > 0) {
|
||||
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
|
||||
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
|
||||
splinePvt->initFromGas(pvdgKeyword, pvtTableIdx, samples);
|
||||
props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
|
||||
} else {
|
||||
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable));
|
||||
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
|
||||
deadPvt->initFromGas(pvdgKeyword, pvtTableIdx);
|
||||
props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
|
||||
}
|
||||
} else if (deck->hasKeyword("PVTG")) {
|
||||
Opm::PvtgTable pvtgTable(deck->getKeyword("PVTG"), /*tableIdx=*/0);
|
||||
|
||||
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable));
|
||||
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG"), pvtTableIdx));
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
|
||||
}
|
||||
@ -102,7 +119,7 @@ std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr dec
|
||||
}
|
||||
|
||||
void testmu(const double reltol, int n, int np, std::vector<double> p, std::vector<double> r,std::vector<double> z,
|
||||
std::vector<std::shared_ptr<SinglePvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
|
||||
std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
|
||||
{
|
||||
std::vector<double> mu(n);
|
||||
std::vector<double> dmudp(n);
|
||||
@ -138,7 +155,7 @@ void testmu(const double reltol, int n, int np, std::vector<double> p, std::vect
|
||||
}
|
||||
|
||||
void testb(const double reltol, int n, int np, std::vector<double> p, std::vector<double> r,std::vector<double> z,
|
||||
std::vector<std::shared_ptr<SinglePvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
|
||||
std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
|
||||
{
|
||||
// test b
|
||||
std::vector<double> b(n);
|
||||
@ -180,7 +197,7 @@ void testb(const double reltol, int n, int np, std::vector<double> p, std::vecto
|
||||
}
|
||||
}
|
||||
|
||||
void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
|
||||
void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<PvtInterface> > props_){
|
||||
// test bublepoint pressure
|
||||
std::vector<double> rs(n);
|
||||
std::vector<double> drsdp(n);
|
||||
@ -202,7 +219,7 @@ void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<
|
||||
}
|
||||
}
|
||||
|
||||
void testrvSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
|
||||
void testrvSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<PvtInterface> > props_){
|
||||
// test rv saturated
|
||||
std::vector<double> rv(n);
|
||||
std::vector<double> drvdp(n);
|
||||
@ -234,9 +251,11 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
|
||||
Opm::ParserPtr parser(new Opm::Parser());
|
||||
Opm::DeckConstPtr deck(parser->parseFile(filename));
|
||||
|
||||
|
||||
|
||||
// setup pvt interface
|
||||
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
|
||||
std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(deck,phase_usage_);
|
||||
std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck,phase_usage_);
|
||||
|
||||
|
||||
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
|
||||
@ -311,7 +330,7 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
|
||||
|
||||
// setup pvt interface
|
||||
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
|
||||
std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(deck,phase_usage_);
|
||||
std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck,phase_usage_);
|
||||
|
||||
|
||||
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
|
||||
|
||||
@ -9,6 +9,19 @@ GAS
|
||||
|
||||
FIELD
|
||||
|
||||
-- tests for the PVT functions need a grid because the OPM-API for the
|
||||
-- PVT functions assumes the presence of compressed grid cells,
|
||||
-- i.e. we need to be able to map from compressed to logical cartesian
|
||||
-- cell indices to find out the PVT region of a cell
|
||||
DIMENS
|
||||
3 3 3 /
|
||||
DXV
|
||||
1 2 3 /
|
||||
DYV
|
||||
1 2 3 /
|
||||
DZV
|
||||
1 2 3 /
|
||||
|
||||
-- PVT PROPERTIES OF DRY GAS (NO VAPOURISED OIL)
|
||||
-- FROM SPE3 Blackoil Kleppe
|
||||
--
|
||||
|
||||
Loading…
Reference in New Issue
Block a user