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

incorperate the review comments/decisions for multi-region PVT

Browse Source 
the largest change is that all classes below opm/core/props/pvt take
the PVT region index as an argument, the higher-level ones (i.e.,
BlackoilProps*) take cell indices.
This commit is contained in:
Andreas Lauser
2014-06-04 11:11:35 +02:00
parent fbd8d42e8d
commit cbb5910a93
27 changed files with 417 additions and 248 deletions
Download Patch File Download Diff File Expand all files Collapse all files
opm/core/pressure/CompressibleTpfa.cpp
+2 -2
View File
@@ -238,7 +238,7 @@ namespace Opm
const int cell = wells_->well_cells[j];
const double cell_depth = grid_.cell_centroids[dim * cell + dim - 1];
props_.matrix(1, &state.pressure()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0);
props_.density(1, &A[0], &rho[0]);
props_.density(1, &A[0], &cell, &rho[0]);
for (int phase = 0; phase < np; ++phase) {
const double s_phase = state.saturation()[np*cell + phase];
wellperf_wdp_[j] += s_phase*rho[phase]*grav*(cell_depth - ref_depth);
@@ -380,7 +380,7 @@ namespace Opm
// Gravity contribution, gravcontrib = rho*(face_z - cell_z) [per phase].
if (grav != 0.0) {
const double depth_diff = face_depth - grid_.cell_centroids[c[j]*dim + dim - 1];
props_.density(1, &cell_A_[np*np*c[j]], &gravcontrib[j][0]);
props_.density(1, &cell_A_[np*np*c[j]], &c[j], &gravcontrib[j][0]);
for (int p = 0; p < np; ++p) {
gravcontrib[j][p] *= depth_diff*grav;
}
opm/core/props/BlackoilPropertiesBasic.cpp
+3 -1
View File
@@ -157,9 +157,11 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
void BlackoilPropertiesBasic::density(const int n,
const double* A,
const int* /*cells*/,
double* rho) const
{
const int np = numPhases();
@@ -177,7 +179,7 @@ namespace Opm
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
const double* BlackoilPropertiesBasic::surfaceDensity() const
const double* BlackoilPropertiesBasic::surfaceDensity(int /*cellIdx*/) const
{
return pvt_.surfaceDensities();
}
opm/core/props/BlackoilPropertiesBasic.hpp
+9 -1
View File
@@ -59,6 +59,11 @@ namespace Opm
/// \return N, the number of cells.
virtual int numCells() const;
/// Return an array containing the PVT table index for each
/// grid cell
virtual const int* cellPvtRegionIndex() const
{ return NULL; }
/// \return Array of N porosity values.
virtual const double* porosity() const;
@@ -114,14 +119,17 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
virtual void density(const int n,
const double* A,
const int* cells,
double* rho) const;
/// Densities of stock components at surface conditions.
/// \param[in] cellIdx The index of the cell for which the surface density ought to be calculated
/// \return Array of P density values.
virtual const double* surfaceDensity() const;
virtual const double* surfaceDensity(int cellIdx = 0) const;
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
opm/core/props/BlackoilPropertiesFromDeck.cpp
+26 -10
View File
@@ -96,14 +96,19 @@ namespace Opm
void BlackoilPropertiesFromDeck::viscosity(const int n,
const double* p,
const double* z,
const int* /*cells*/,
const int* cells,
double* mu,
double* dmudp) const
{
if (dmudp) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::viscosity() -- derivatives of viscosity not yet implemented.");
} else {
pvt_.mu(n, p, z, mu);
const int *cellPvtTableIdx = cellPvtRegionIndex();
std::vector<int> pvtTableIdx(n);
for (int i = 0; i < n; ++ i)
pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
pvt_.mu(n, &pvtTableIdx[0], p, z, mu);
}
}
@@ -120,21 +125,27 @@ namespace Opm
void BlackoilPropertiesFromDeck::matrix(const int n,
const double* p,
const double* z,
const int* /*cells*/,
const int* cells,
double* A,
double* dAdp) const
{
const int np = numPhases();
const int *cellPvtTableIdx = cellPvtRegionIndex();
std::vector<int> pvtTableIdx(n);
for (int i = 0; i < n; ++ i)
pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
B_.resize(n*np);
R_.resize(n*np);
if (dAdp) {
dB_.resize(n*np);
dR_.resize(n*np);
pvt_.dBdp(n, p, z, &B_[0], &dB_[0]);
pvt_.dRdp(n, p, z, &R_[0], &dR_[0]);
pvt_.dBdp(n, &pvtTableIdx[0], p, z, &B_[0], &dB_[0]);
pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]);
} else {
pvt_.B(n, p, z, &B_[0]);
pvt_.R(n, p, z, &R_[0]);
pvt_.B(n, &pvtTableIdx[0], p, z, &B_[0]);
pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]);
}
const int* phase_pos = pvt_.phasePosition();
bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] &&
@@ -207,15 +218,19 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
void BlackoilPropertiesFromDeck::density(const int n,
const double* A,
const int* cells,
double* rho) const
{
const int np = numPhases();
const double* sdens = pvt_.surfaceDensities();
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int cellIdx = cells?cells[i]:i;
int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx);
const double* sdens = pvt_.surfaceDensities(pvtRegionIdx);
for (int phase = 0; phase < np; ++phase) {
rho[np*i + phase] = 0.0;
for (int comp = 0; comp < np; ++comp) {
@@ -227,9 +242,10 @@ namespace Opm
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
const double* BlackoilPropertiesFromDeck::surfaceDensity() const
const double* BlackoilPropertiesFromDeck::surfaceDensity(int cellIdx) const
{
return pvt_.surfaceDensities();
int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx);
return pvt_.surfaceDensities(pvtRegionIdx);
}
/// \param[in] n Number of data points.
opm/core/props/BlackoilPropertiesFromDeck.hpp
+16 -1
View File
@@ -97,6 +97,11 @@ namespace Opm
/// \return N, the number of cells.
virtual int numCells() const;
/// Return an array containing the PVT table index for each
/// grid cell
virtual const int* cellPvtRegionIndex() const
{ return &cellPvtRegionIdx_[0]; }
/// \return Array of N porosity values.
virtual const double* porosity() const;
@@ -152,14 +157,16 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
virtual void density(const int n,
const double* A,
const int* cells,
double* rho) const;
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
virtual const double* surfaceDensity() const;
virtual const double* surfaceDensity(int cellIdx = 0) const;
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
@@ -206,6 +213,13 @@ namespace Opm
double* smax) const;
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
template<class CentroidIterator>
void init(Opm::DeckConstPtr deck,
int number_of_cells,
@@ -224,6 +238,7 @@ namespace Opm
const parameter::ParameterGroup& param,
bool init_rock);
RockFromDeck rock_;
std::vector<int> cellPvtRegionIdx_;
BlackoilPvtProperties pvt_;
std::unique_ptr<SaturationPropsInterface> satprops_;
mutable std::vector<double> B_;
opm/core/props/BlackoilPropertiesFromDeck_impl.hpp
+9 -2
View File
@@ -42,10 +42,14 @@ namespace Opm
int dimension,
bool init_rock)
{
// retrieve the cell specific PVT table index from the deck
// and using the grid...
extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell);
if (init_rock){
rock_.init(deck, number_of_cells, global_cell, cart_dims);
}
pvt_.init(deck, /*numSamples=*/0, number_of_cells, global_cell);
pvt_.init(deck, /*numSamples=*/0);
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
@@ -68,13 +72,16 @@ namespace Opm
const parameter::ParameterGroup& param,
bool init_rock)
{
// retrieve the cell specific PVT table index from the deck
// and using the grid...
extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell);
if(init_rock){
rock_.init(deck, number_of_cells, global_cell, cart_dims);
}
const int pvt_samples = param.getDefault("pvt_tab_size", 200);
pvt_.init(deck, pvt_samples, number_of_cells, global_cell);
pvt_.init(deck, pvt_samples);
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 200);
opm/core/props/BlackoilPropertiesInterface.hpp
+7 -1
View File
@@ -47,6 +47,10 @@ namespace Opm
/// \return N, the number of cells.
virtual int numCells() const = 0;
/// Return an array containing the PVT table index for each
/// grid cell
virtual const int* cellPvtRegionIndex() const = 0;
/// \return Array of N porosity values.
virtual const double* porosity() const = 0;
@@ -102,14 +106,16 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
virtual void density(const int n,
const double* A,
const int* cells,
double* rho) const = 0;
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
virtual const double* surfaceDensity() const = 0;
virtual const double* surfaceDensity(int regionIdx = 0) const = 0;
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
opm/core/props/pvt/BlackoilPvtProperties.cpp
+19 -24
View File
@@ -43,9 +43,7 @@ namespace Opm
}
void BlackoilPvtProperties::init(Opm::DeckConstPtr deck,
int numSamples,
int numCompressedCells,
const int *compressedToCartesianCellIdx)
int numSamples)
{
phase_usage_ = phaseUsageFromDeck(deck);
@@ -69,14 +67,6 @@ namespace Opm
}
}
// 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);
@@ -85,7 +75,7 @@ namespace Opm
// 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_);
pvtw->initFromWater(deck->getKeyword("PVTW"));
props_[phase_usage_.phase_pos[Aqua]] = pvtw;
}
@@ -97,18 +87,18 @@ namespace Opm
Opm::DeckKeywordConstPtr pvdoKeyword = deck->getKeyword("PVDO");
if (numSamples > 0) {
auto splinePvt = std::shared_ptr<PvtDeadSpline>(new PvtDeadSpline);
splinePvt->initFromOil(pvdoKeyword, pvtTableIdx_, numSamples);
splinePvt->initFromOil(pvdoKeyword, numSamples);
props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
} else {
auto deadPvt = std::shared_ptr<PvtDead>(new PvtDead);
deadPvt->initFromOil(pvdoKeyword, pvtTableIdx_);
deadPvt->initFromOil(pvdoKeyword);
props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
}
} else if (deck->hasKeyword("PVTO")) {
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"), pvtTableIdx_));
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO")));
} else if (deck->hasKeyword("PVCDO")) {
std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
pvcdo->initFromOil(deck->getKeyword("PVCDO"), pvtTableIdx_);
pvcdo->initFromOil(deck->getKeyword("PVCDO"));
props_[phase_usage_.phase_pos[Liquid]] = pvcdo;
} else {
@@ -123,17 +113,17 @@ namespace Opm
if (numSamples > 0) {
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromGas(pvdgKeyword, pvtTableIdx_, numSamples);
splinePvt->initFromGas(pvdgKeyword, numSamples);
props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
} else {
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromGas(pvdgKeyword, pvtTableIdx_);
deadPvt->initFromGas(pvdgKeyword);
props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
}
} else if (deck->hasKeyword("PVTG")) {
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG"), pvtTableIdx_));
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG")));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
}
@@ -168,13 +158,14 @@ namespace Opm
void BlackoilPvtProperties::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const
{
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->mu(n, p, z, &data1_[0]);
props_[phase]->mu(n, pvtTableIdx, p, z, &data1_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[phase_usage_.num_phases*i + phase] = data1_[i];
@@ -183,13 +174,14 @@ namespace Opm
}
void BlackoilPvtProperties::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B) const
{
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->B(n, p, z, &data1_[0]);
props_[phase]->B(n, pvtTableIdx, p, z, &data1_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i];
@@ -198,6 +190,7 @@ namespace Opm
}
void BlackoilPvtProperties::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@@ -206,7 +199,7 @@ namespace Opm
data1_.resize(n);
data2_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->dBdp(n, p, z, &data1_[0], &data2_[0]);
props_[phase]->dBdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i];
@@ -217,13 +210,14 @@ namespace Opm
void BlackoilPvtProperties::R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const
{
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->R(n, p, z, &data1_[0]);
props_[phase]->R(n, pvtTableIdx, p, z, &data1_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[phase_usage_.num_phases*i + phase] = data1_[i];
@@ -232,6 +226,7 @@ namespace Opm
}
void BlackoilPvtProperties::dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
@@ -240,7 +235,7 @@ namespace Opm
data1_.resize(n);
data2_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->dRdp(n, p, z, &data1_[0], &data2_[0]);
props_[phase]->dRdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[phase_usage_.num_phases*i + phase] = data1_[i];
opm/core/props/pvt/BlackoilPvtProperties.hpp
+15 -23
View File
@@ -33,10 +33,12 @@ namespace Opm
{
/// Class collecting the pvt properties for all active phases.
/// For all the methods, the following apply: p and z
/// are expected to be of size n and n*num_phases, respectively.
/// Output arrays shall be of size n*num_phases, and must be valid
/// before calling the method.
/// For all the methods, the following apply:
/// - p and z are expected to be of size n and n*num_phases, respectively.
/// - pvtTableIdx specifies the PVT table to be used for each data
/// point and is thus expected to be an array of size n
/// - 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 PvtInterface is that this collects all phases' properties,
/// and therefore the output arrays are of size n*num_phases as opposed
@@ -49,16 +51,9 @@ namespace Opm
/// 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,
int numCompressedCells,
const int *compressedToCartesianCellIdx);
int samples);
/// \return Object describing the active phases.
PhaseUsage phaseUsage() const;
@@ -82,18 +77,21 @@ namespace Opm
/// Viscosity as a function of p and z.
void mu(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const;
/// Formation volume factor as a function of p and z.
void B(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
void dBdp(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@@ -101,12 +99,14 @@ namespace Opm
/// Solution factor as a function of p and z.
void R(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
void dRdp(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_R,
@@ -119,19 +119,11 @@ namespace Opm
PhaseUsage phase_usage_;
// The PVT properties. One object per active fluid phase.
// The PVT properties. We need to store one object per PVT
// region per active fluid phase.
std::vector<std::shared_ptr<PvtInterface> > 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_;
std::vector<std::array<double, MaxNumPhases> > densities_;
mutable std::vector<double> data1_;
mutable std::vector<double> data2_;
};
opm/core/props/pvt/PvtConstCompr.hpp
+40 -24
View File
@@ -34,22 +34,24 @@ 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.
/// 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). Also, since this class supports
/// multiple PVT regions, the concrete table to be used for each
/// data point needs to be specified via the pvtTableIdx argument
/// of the respective method. For all the virtual methods, the
/// following apply: pvtTableIdx, p, r and z are expected to be of
/// size n, 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)
void initFromWater(Opm::DeckKeywordConstPtr pvtwKeyword)
{
pvtTableIdx_ = pvtTableIdx;
int numRegions = pvtwKeyword->size();
ref_press_.resize(numRegions);
@@ -69,10 +71,8 @@ namespace Opm
}
}
void initFromOil(Opm::DeckKeywordConstPtr pvcdoKeyword, const std::vector<int> &pvtTableIdx)
void initFromOil(Opm::DeckKeywordConstPtr pvcdoKeyword)
{
pvtTableIdx_ = pvtTableIdx;
int numRegions = pvcdoKeyword->size();
ref_press_.resize(numRegions);
@@ -92,7 +92,11 @@ namespace Opm
}
}
PvtConstCompr(double visc)
/*!
* \brief Create a PVT object with a given viscosity that
* assumes all fluid phases to be incompressible.
*/
explicit PvtConstCompr(double visc)
: ref_press_(1, 0.0),
ref_B_(1, 1.0),
comp_(1, 0.0),
@@ -106,6 +110,7 @@ namespace Opm
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*z*/,
double* output_mu) const
@@ -113,13 +118,14 @@ namespace Opm
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(i);
int tableIdx = getTableIndex_(pvtRegionIdx, 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 int* pvtRegionIdx,
const double* p,
const double* /*r*/,
double* output_mu,
@@ -129,7 +135,7 @@ namespace Opm
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(i);
int tableIdx = getTableIndex_(pvtRegionIdx, 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;
@@ -139,6 +145,7 @@ namespace Opm
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
@@ -149,7 +156,7 @@ namespace Opm
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(i);
int tableIdx = getTableIndex_(pvtRegionIdx, 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;
@@ -159,6 +166,7 @@ namespace Opm
}
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*z*/,
double* output_B) const
@@ -166,13 +174,14 @@ namespace Opm
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(i);
int tableIdx = getTableIndex_(pvtRegionIdx, 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 int* pvtRegionIdx,
const double* p,
const double* /*z*/,
double* output_B,
@@ -180,7 +189,7 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(i);
int tableIdx = getTableIndex_(pvtRegionIdx, 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;
@@ -189,6 +198,7 @@ namespace Opm
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*r*/,
double* output_b,
@@ -198,7 +208,7 @@ namespace Opm
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(i);
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
@@ -211,6 +221,7 @@ namespace Opm
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
@@ -221,7 +232,7 @@ namespace Opm
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(i);
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
@@ -233,6 +244,7 @@ namespace Opm
}
virtual void rsSat(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
@@ -242,6 +254,7 @@ namespace Opm
}
virtual void rvSat(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
@@ -251,6 +264,7 @@ namespace Opm
}
virtual void R(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
@@ -259,6 +273,7 @@ namespace Opm
}
virtual void dRdp(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R,
@@ -269,14 +284,15 @@ namespace Opm
}
private:
int getTableIndex_(int cellIdx) const
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (pvtTableIdx_.empty())
if (!pvtTableIdx)
return 0;
return pvtTableIdx_[cellIdx];
return pvtTableIdx[cellIdx];
}
std::vector<int> pvtTableIdx_;
// The PVT properties. We need to store one value per PVT
// region.
std::vector<double> ref_press_;
std::vector<double> ref_B_;
std::vector<double> comp_;
opm/core/props/pvt/PvtDead.cpp
+21 -16
View File
@@ -34,11 +34,8 @@ namespace Opm
// Member functions
//-------------------------------------------------------------------------
/// Constructor
void PvtDead::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
const std::vector<int> &pvtTableIdx)
void PvtDead::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword)
{
pvtTableIdx_ = pvtTableIdx;
int numRegions = Opm::PvdoTable::numTables(pvdoKeyword);
// resize the attributes of the object
@@ -63,11 +60,8 @@ namespace Opm
}
void PvtDead::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
const std::vector<int> &pvtTableIdx)
void PvtDead::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword)
{
pvtTableIdx_ = pvtTableIdx;
int numRegions = Opm::PvdgTable::numTables(pvdgKeyword);
// resize the attributes of the object
@@ -100,18 +94,20 @@ namespace Opm
void PvtDead::mu(const int n,
const int* pvtTableIdx,
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];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
}
}
void PvtDead::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
double* output_mu,
@@ -120,7 +116,7 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
}
@@ -129,6 +125,7 @@ namespace Opm
}
void PvtDead::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
@@ -138,7 +135,7 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
}
@@ -147,6 +144,7 @@ namespace Opm
}
void PvtDead::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*z*/,
double* output_B) const
@@ -154,27 +152,29 @@ namespace Opm
// #pragma omp parallel for
// B = 1/b
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_B[i] = 1.0/b_[regionIdx](p[i]);
}
}
void PvtDead::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
B(n, p, 0, output_B);
B(n, pvtTableIdx, p, 0, output_B);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
double Bg = output_B[i];
output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]);
}
}
void PvtDead::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
double* output_b,
@@ -184,7 +184,7 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
@@ -195,6 +195,7 @@ namespace Opm
}
void PvtDead::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
@@ -205,7 +206,7 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
@@ -216,6 +217,7 @@ namespace Opm
}
void PvtDead::rsSat(const int n,
const int* pvtTableIdx,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
@@ -225,6 +227,7 @@ namespace Opm
}
void PvtDead::rvSat(const int n,
const int* pvtTableIdx,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
@@ -234,6 +237,7 @@ namespace Opm
}
void PvtDead::R(const int n,
const int* pvtTableIdx,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
@@ -242,6 +246,7 @@ namespace Opm
}
void PvtDead::dRdp(const int n,
const int* pvtTableIdx,
const double* /*p*/,
const double* /*z*/,
double* output_R,
opm/core/props/pvt/PvtDead.hpp
+32 -12
View File
@@ -33,25 +33,28 @@ namespace Opm
{
/// Class for immiscible dead oil and dry gas.
/// 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.
/// 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). Also, since this class supports
/// multiple PVT regions, the concrete table to be used for each
/// data point needs to be specified via the pvtTableIdx argument
/// of the respective method. For all the virtual methods, the
/// following apply: pvtTableIdx, p, r and z are expected to be of
/// size n, 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 PvtDead : public PvtInterface
{
public:
PvtDead() {};
void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
const std::vector<int> &pvtTableIdx);
void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
const std::vector<int> &pvtTableIdx);
void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword);
void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword);
virtual ~PvtDead();
/// Viscosity as a function of p and z.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const;
@@ -59,6 +62,7 @@ namespace Opm
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_mu,
@@ -68,6 +72,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -77,12 +82,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@@ -91,6 +98,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_b,
@@ -100,6 +108,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -110,31 +119,42 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
private:
// PVT properties of dry gas or dead oil
std::vector<int> pvtTableIdx_;
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
// PVT properties of dry gas or dead oil. We need to store one
// table per PVT region.
std::vector<NonuniformTableLinear<double> > b_;
std::vector<NonuniformTableLinear<double> > viscosity_;
};
opm/core/props/pvt/PvtDeadSpline.cpp
+19 -14
View File
@@ -40,11 +40,8 @@ namespace Opm
{}
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
@@ -72,11 +69,8 @@ namespace Opm
}
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
@@ -111,18 +105,20 @@ namespace Opm
void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx,
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];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
}
}
void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
double* output_mu,
@@ -131,7 +127,7 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
}
@@ -139,6 +135,7 @@ namespace Opm
}
void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
@@ -149,7 +146,7 @@ namespace Opm
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
}
@@ -157,33 +154,36 @@ namespace Opm
}
void PvtDeadSpline::B(const int n,
const int* pvtTableIdx,
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];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_B[i] = 1.0/b_[regionIdx](p[i]);
}
}
void PvtDeadSpline::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
B(n, p, 0, output_B);
B(n, pvtTableIdx, p, 0, output_B);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
double Bg = output_B[i];
output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]);
}
}
void PvtDeadSpline::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
double* output_b,
@@ -192,7 +192,7 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
}
@@ -200,6 +200,7 @@ namespace Opm
}
void PvtDeadSpline::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*r*/,
const PhasePresence* /*cond*/,
@@ -209,7 +210,7 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = pvtTableIdx_[i];
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
}
@@ -217,6 +218,7 @@ namespace Opm
}
void PvtDeadSpline::rsSat(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
@@ -226,6 +228,7 @@ namespace Opm
}
void PvtDeadSpline::rvSat(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
@@ -235,6 +238,7 @@ namespace Opm
}
void PvtDeadSpline::R(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
@@ -243,6 +247,7 @@ namespace Opm
}
void PvtDeadSpline::dRdp(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R,
opm/core/props/pvt/PvtDeadSpline.hpp
+20 -4
View File
@@ -44,16 +44,15 @@ namespace Opm
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,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const;
@@ -61,6 +60,7 @@ namespace Opm
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_mu,
@@ -70,6 +70,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -79,12 +80,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@@ -93,6 +96,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_b,
@@ -102,6 +106,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -111,31 +116,42 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
private:
// PVT properties of dry gas or dead oil
std::vector<int> pvtTableIdx_;
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
// PVT properties of dry gas or dead oil. We need to store one
// table per PVT region.
std::vector<UniformTableLinear<double> > b_;
std::vector<UniformTableLinear<double> > viscosity_;
};
opm/core/props/pvt/PvtInterface.hpp
+31 -12
View File
@@ -44,13 +44,19 @@ namespace Opm
/// 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.
/// For all the virtual methods, the following apply:
/// - pvtRegionIdx is an array of size n and represents the
/// index of the PVT table which should be used to calculate
/// the output. NULL can also be passed and is interpreted
/// such that the first table should be used for the output
/// - 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.
/// Viscosity as a function of p and z.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_mu) const = 0;
@@ -58,6 +64,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
double* output_mu,
@@ -67,6 +74,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -76,12 +84,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B) const = 0;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B,
@@ -90,6 +100,7 @@ namespace Opm
/// The inverse of the volume factor b = 1 / B as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
double* output_b,
@@ -99,6 +110,7 @@ namespace Opm
/// The inverse of the volume factor b = 1 / B as a function of p and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -108,12 +120,14 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const = 0;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const = 0;
@@ -121,12 +135,14 @@ namespace Opm
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const = 0;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
@@ -140,18 +156,21 @@ namespace Opm
/*!
* \brief Helper function to create an array containing the (C-Style)
* PVT table index for each compressed cell.
* PVT table index for each compressed cell from an Eclipse deck.
*
* 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.
* This function assumes that the degrees of freedom where PVT
* properties need to be calculated are grid cells. 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 while the classes which use the PvtInterface
* implementations usually use compressed cells.
*/
void extractPvtTableIndex(std::vector<int> &pvtTableIdx,
void extractPvtTableIndex(std::vector<int>& pvtTableIdx,
Opm::DeckConstPtr deck,
size_t numCompressed,
const int *compressedToCartesianIdx);
const int* compressedToCartesianIdx);
} // namespace Opm
opm/core/props/pvt/PvtLiveGas.cpp
+26 -16
View File
@@ -46,11 +46,8 @@ namespace Opm
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
PvtLiveGas::PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword,
const std::vector<int> &pvtTableIdx)
PvtLiveGas::PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword)
{
pvtTableIdx_ = pvtTableIdx;
int numTables = Opm::PvtgTable::numTables(pvtgKeyword);
saturated_gas_table_.resize(numTables);
undersat_gas_tables_.resize(numTables);
@@ -93,18 +90,20 @@ namespace Opm
void PvtLiveGas::mu(const int n,
const int* pvtRegionIdx,
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);
output_mu[i] = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 2, false);
}
}
/// Viscosity and its derivatives as a function of p and r.
void PvtLiveGas::mu(const int /*n*/,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*r*/,
double* /*output_mu*/,
@@ -116,6 +115,7 @@ namespace Opm
/// Viscosity and its derivatives as a function of p and r.
void PvtLiveGas::mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -125,9 +125,10 @@ namespace Opm
{
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);
int tableIdx = getTableIndex_(pvtRegionIdx, i);
output_mu[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 0);
output_dmudp[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 1);
output_dmudr[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 2);
}
}
@@ -135,13 +136,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
void PvtLiveGas::B(const int n,
const int* pvtRegionIdx,
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]);
output_B[i] = evalB(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i));
}
}
@@ -149,6 +151,7 @@ namespace Opm
/// Formation volume factor and p-derivative as functions of p and z.
void PvtLiveGas::dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B,
@@ -156,12 +159,13 @@ namespace Opm
{
// #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]);
evalBDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, 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 int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*r*/,
double* /*output_b*/,
@@ -174,6 +178,7 @@ namespace Opm
/// 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 int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -186,21 +191,23 @@ namespace Opm
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);
int tableIdx = getTableIndex_(pvtRegionIdx, i);
output_b[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 0);
output_dbdp[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 1);
output_dbdr[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 2);
}
}
/// Gas resolution and its derivatives at bublepoint as a function of p.
void PvtLiveGas::rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const
{
for (int i = 0; i < n; ++i) {
int pvtTableIdx = pvtTableIdx_[i];
int pvtTableIdx = getTableIndex_(pvtRegionIdx, 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],
@@ -210,6 +217,7 @@ namespace Opm
}
void PvtLiveGas::rsSat(const int n,
const int* pvtRegionIdx,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
@@ -220,13 +228,14 @@ namespace Opm
/// Solution factor as a function of p and z.
void PvtLiveGas::R(const int n,
const int* pvtRegionIdx,
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]);
output_R[i] = evalR(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i));
}
}
@@ -234,6 +243,7 @@ namespace Opm
/// Solution factor and p-derivative as functions of p and z.
void PvtLiveGas::dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
@@ -241,7 +251,7 @@ namespace Opm
{
// #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]);
evalRDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_R[i], output_dRdp[i]);
}
}
opm/core/props/pvt/PvtLiveGas.hpp
+21 -3
View File
@@ -38,11 +38,12 @@ namespace Opm
class PvtLiveGas : public PvtInterface
{
public:
PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword, const std::vector<int> &pvtTableIdx);
PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword);
virtual ~PvtLiveGas();
/// Viscosity as a function of p and z.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_mu) const;
@@ -50,6 +51,7 @@ namespace Opm
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
double* output_mu,
@@ -59,6 +61,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -68,12 +71,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_B,
@@ -82,6 +87,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
double* output_b,
@@ -91,6 +97,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -102,30 +109,41 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
protected:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
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;
@@ -143,8 +161,8 @@ namespace Opm
const int pvtTableIdx,
const int item,
const int deriv = 0) const;
// PVT properties of wet gas (with vaporised oil)
std::vector<int> pvtTableIdx_;
// PVT properties of wet gas (with vaporised oil). We need to
// store one table per PVT region.
std::vector< std::vector<std::vector<double> > > saturated_gas_table_;
std::vector< std::vector<std::vector<std::vector<double> > > > undersat_gas_tables_;
};
opm/core/props/pvt/PvtLiveOil.cpp
+43 -35
View File
@@ -34,11 +34,8 @@ namespace Opm
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
PvtLiveOil::PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword,
const std::vector<int> &pvtTableIdx)
PvtLiveOil::PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword)
{
pvtTableIdx_ = pvtTableIdx;
int numTables = Opm::PvtoTable::numTables(pvtoKeyword);
saturated_oil_table_.resize(numTables);
undersat_oil_tables_.resize(numTables);
@@ -120,20 +117,22 @@ namespace Opm
/// Viscosity as a function of p and z.
void PvtLiveOil::mu(const int n,
const int* pvtTableIdx,
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];
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = miscible_oil(p[i], z + num_phases_*i, pvtTableIdx, 2, false);
output_mu[i] = miscible_oil(p[i], z + num_phases_*i, tableIdx, 2, false);
}
}
/// Viscosity and its derivatives as a function of p and r.
void PvtLiveOil::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_mu,
@@ -142,17 +141,18 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int pvtTableIdx = pvtTableIdx_[i];
int tableIdx = getTableIndex_(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);
output_mu[i] = miscible_oil(p[i], r[i], tableIdx, 2, 0);
output_dmudp[i] = miscible_oil(p[i], r[i], tableIdx, 2, 1);
output_dmudr[i] = miscible_oil(p[i], r[i], tableIdx, 2, 2);
}
}
/// Viscosity and its derivatives as a function of p and r.
void PvtLiveOil::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -162,12 +162,12 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int pvtTableIdx = pvtTableIdx_[i];
int tableIdx = getTableIndex_(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);
output_mu[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 0);
output_dmudp[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 1);
output_dmudr[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 2);
}
}
@@ -175,15 +175,16 @@ namespace Opm
/// Formation volume factor as a function of p and z.
void PvtLiveOil::B(const int n,
const int* pvtTableIdx,
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];
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_B[i] = evalB(pvtTableIdx, p[i], z + num_phases_*i);
output_B[i] = evalB(tableIdx, p[i], z + num_phases_*i);
}
}
@@ -191,6 +192,7 @@ namespace Opm
/// Formation volume factor and p-derivative as functions of p and z.
void PvtLiveOil::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@@ -198,13 +200,14 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int pvtTableIdx = pvtTableIdx_[i];
int tableIdx = getTableIndex_(pvtTableIdx, i);
evalBDeriv(pvtTableIdx, p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
evalBDeriv(tableIdx, p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
}
}
void PvtLiveOil::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_b,
@@ -214,16 +217,17 @@ namespace Opm
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int pvtTableIdx = pvtTableIdx_[i];
int tableIdx = getTableIndex_(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);
output_b[i] = miscible_oil(tableIdx, p[i], r[i], 1, 0);
output_dbdp[i] = miscible_oil(tableIdx, p[i], r[i], 1, 1);
output_dbdr[i] = miscible_oil(tableIdx, p[i], r[i], 1, 2);
}
}
void PvtLiveOil::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -235,32 +239,34 @@ namespace Opm
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
int pvtTableIdx = pvtTableIdx_[i];
int tableIdx = getTableIndex_(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);
output_b[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 0);
output_dbdp[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 1);
output_dbdr[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 2);
}
}
void PvtLiveOil::rsSat(const int n,
const int* pvtTableIdx,
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]);
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_rsSat[i] = linearInterpolation(saturated_oil_table_[tableIdx][0],
saturated_oil_table_[tableIdx][3],p[i]);
output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[tableIdx][0],
saturated_oil_table_[tableIdx][3],p[i]);
}
}
void PvtLiveOil::rvSat(const int n,
const int* pvtTableIdx,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
@@ -271,14 +277,15 @@ namespace Opm
/// Solution factor as a function of p and z.
void PvtLiveOil::R(const int n,
const int* pvtTableIdx,
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);
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_R[i] = evalR(tableIdx, p[i], z + num_phases_*i);
}
}
@@ -286,6 +293,7 @@ namespace Opm
/// Solution factor and p-derivative as functions of p and z.
void PvtLiveOil::dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
@@ -293,8 +301,8 @@ namespace Opm
{
// #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]);
int tableIdx = getTableIndex_(pvtTableIdx, i);
evalRDeriv(tableIdx, p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
}
}
opm/core/props/pvt/PvtLiveOil.hpp
+21 -3
View File
@@ -39,11 +39,12 @@ namespace Opm
class PvtLiveOil : public PvtInterface
{
public:
PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword, const std::vector<int> &pvtTableIdx);
PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword);
virtual ~PvtLiveOil();
/// Viscosity as a function of p and z.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_mu) const;
@@ -51,6 +52,7 @@ namespace Opm
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_mu,
@@ -60,6 +62,7 @@ namespace Opm
/// Viscosity as a function of p and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -69,12 +72,14 @@ namespace Opm
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_B,
@@ -83,6 +88,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
double* output_b,
@@ -92,6 +98,7 @@ namespace Opm
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* r,
const PhasePresence* cond,
@@ -101,30 +108,41 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
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;
@@ -150,8 +168,8 @@ namespace Opm
const int item,
const int deriv = 0) const;
// PVT properties of live oil (with dissolved gas)
std::vector<int> pvtTableIdx_;
// PVT properties of live oil (with dissolved gas). We need to
// store one table per PVT region.
std::vector<std::vector<std::vector<double> > > saturated_oil_table_;
std::vector<std::vector<std::vector<std::vector<double> > > > undersat_oil_tables_;
};
opm/core/props/pvt/PvtPropertiesBasic.hpp
+2
View File
@@ -95,6 +95,8 @@ namespace Opm
double* output_dRdp) const;
private:
// The PVT properties. We need to store one value per PVT
// region.
std::vector<double> density_;
std::vector<double> viscosity_;
std::vector<double> formation_volume_factor_;
opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
+2 -2
View File
@@ -33,9 +33,9 @@ namespace Opm
{
}
void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck )
void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck)
{
// If we need multiple regions, this class and the Pvt* classes must change.
// So far, this class only supports a single PVT region. TODO?
int region_number = 0;
PhaseUsage phase_usage = phaseUsageFromDeck(deck);
opm/core/simulator/EquilibrationHelpers.hpp
+1 -1
View File
@@ -139,7 +139,7 @@ namespace Opm
props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
std::vector<double> rho(np, 0.0);
props_.density(1, &A[0], &rho[0]);
props_.density(1, &A[0], &c_[0], &rho[0]);
return rho;
}
opm/core/simulator/initStateEquil.hpp
-2
View File
@@ -369,8 +369,6 @@ namespace Opm
const UnstructuredGrid& G ,
const double grav)
{
typedef Miscibility::NoMixing NoMix;
for (typename RMap::RegionId
r = 0, nr = reg.numRegions();
r < nr; ++r)
opm/core/simulator/initState_impl.hpp
+1 -1
View File
@@ -184,7 +184,7 @@ namespace Opm
double A[4] = { 0.0 };
props_.matrix(1, &pressure, surfvol[phase], cells, A, 0);
double rho[2] = { 0.0 };
props_.density(1, A, rho);
props_.density(1, A, cells, rho);
return rho[phase];
}
};
opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp
+1 -1
View File
@@ -418,7 +418,7 @@ namespace Opm
assert(np == 2);
const int dim = grid_.dimensions;
density_.resize(nc*np);
props_.density(grid_.number_of_cells, &A_[0], &density_[0]);
props_.density(grid_.number_of_cells, &A_[0], /*cellIndices=*/NULL, &density_[0]);
std::fill(gravflux_.begin(), gravflux_.end(), 0.0);
for (int f = 0; f < nf; ++f) {
const int* c = &grid_.face_cells[2*f];
opm/core/utility/miscUtilitiesBlackoil.cpp
+2 -1
View File
@@ -285,7 +285,8 @@ namespace Opm
/// @brief Computes saturation from surface volume
void computeSaturation(const BlackoilPropertiesInterface& props,
BlackoilState& state){
BlackoilState& state)
{
const int np = props.numPhases();
const int nc = props.numCells();
tests/test_blackoilfluid.cpp
+28 -36
View File
@@ -40,20 +40,10 @@ using namespace std;
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<PvtInterface> > props_;
// Set the properties.
props_.resize(phase_usage_.num_phases);
@@ -62,7 +52,7 @@ std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, Pha
if (phase_usage_.phase_used[Aqua]) {
if (deck->hasKeyword("PVTW")) {
std::shared_ptr<PvtConstCompr> pvtw(new PvtConstCompr);
pvtw->initFromWater(deck->getKeyword("PVTW"), pvtTableIdx);
pvtw->initFromWater(deck->getKeyword("PVTW"));
props_[phase_usage_.phase_pos[Aqua]] = pvtw;
} else {
@@ -77,18 +67,18 @@ std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, Pha
Opm::DeckKeywordConstPtr pvdoKeyword(deck->getKeyword("PVDO"));
if (samples > 0) {
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromOil(pvdoKeyword, pvtTableIdx, samples);
splinePvt->initFromOil(pvdoKeyword, samples);
props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
} else {
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromOil(pvdoKeyword, pvtTableIdx);
deadPvt->initFromOil(pvdoKeyword);
props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
}
} else if (deck->hasKeyword("PVTO")) {
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO"), pvtTableIdx));
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO")));
} else if (deck->hasKeyword("PVCDO")) {
std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
pvcdo->initFromOil(deck->getKeyword("PVCDO"), pvtTableIdx);
pvcdo->initFromOil(deck->getKeyword("PVCDO"));
props_[phase_usage_.phase_pos[Liquid]] = pvcdo;
} else {
@@ -101,15 +91,15 @@ std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, Pha
Opm::DeckKeywordConstPtr pvdgKeyword(deck->getKeyword("PVDG"));
if (samples > 0) {
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromGas(pvdgKeyword, pvtTableIdx, samples);
splinePvt->initFromGas(pvdgKeyword, samples);
props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
} else {
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromGas(pvdgKeyword, pvtTableIdx);
deadPvt->initFromGas(pvdgKeyword);
props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
}
} else if (deck->hasKeyword("PVTG")) {
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG"), pvtTableIdx));
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG")));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
}
@@ -118,7 +108,7 @@ std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, Pha
return props_;
}
void testmu(const double reltol, int n, int np, std::vector<double> p, std::vector<double> r,std::vector<double> z,
void testmu(const double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{
std::vector<double> mu(n);
@@ -132,8 +122,8 @@ void testmu(const double reltol, int n, int np, std::vector<double> p, std::vect
// test mu
for (int phase = 0; phase < np; ++phase) {
props_[phase]->mu(n, &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
props_[phase]->mu(n, &p[0], &z[0], &mu[0]);
props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &z[0], &mu[0]);
dmudp_diff = (mu_new[1]-mu_new[0])/(p[1]-p[0]);
dmudr_diff = (mu_new[2]-mu_new[0])/(r[2]-r[0]);
dmudp_diff_u = (mu_new[4]-mu_new[3])/(p[4]-p[3]);
@@ -154,7 +144,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,
void testb(const double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{
// test b
@@ -171,8 +161,8 @@ void testb(const double reltol, int n, int np, std::vector<double> p, std::vecto
double dbdr_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase]->b(n, &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
props_[phase]->dBdp(n, &p[0], &z[0], &B[0], &dBdp[0]);
props_[phase]->b(n, &pvtTableIdx[0], &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
props_[phase]->dBdp(n, &pvtTableIdx[0], &p[0], &z[0], &B[0], &dBdp[0]);
dbdp_diff = (b[1]-b[0])/(p[1]-p[0]);
dbdr_diff = (b[2]-b[0])/(r[2]-r[0]);
dbdp_diff_u = (b[4]-b[3])/(p[4]-p[3]);
@@ -197,7 +187,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<PvtInterface> > props_){
void testrsSat(double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<std::shared_ptr<PvtInterface> > props_){
// test bublepoint pressure
std::vector<double> rs(n);
std::vector<double> drsdp(n);
@@ -205,7 +195,7 @@ void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<
double drsdp_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase] ->rsSat(n, &p[0], &rs[0], &drsdp[0]);
props_[phase] ->rsSat(n, &pvtTableIdx[0], &p[0], &rs[0], &drsdp[0]);
drsdp_diff = (rs[1]-rs[0])/(p[1]-p[0]);
drsdp_diff_u = (rs[4]-rs[3])/(p[4]-p[3]);
@@ -219,7 +209,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<PvtInterface> > props_){
void testrvSat(double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<std::shared_ptr<PvtInterface> > props_){
// test rv saturated
std::vector<double> rv(n);
std::vector<double> drvdp(n);
@@ -227,7 +217,7 @@ void testrvSat(double reltol, int n, int np, std::vector<double> p, std::vector<
double drvdp_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase] ->rvSat(n, &p[0], &rv[0], &drvdp[0]);
props_[phase] ->rvSat(n, &pvtTableIdx[0], &p[0], &rv[0], &drvdp[0]);
drvdp_diff = (rv[1]-rv[0])/(p[1]-p[0]);
drvdp_diff_u = (rv[4]-rv[3])/(p[4]-p[3]);
@@ -262,6 +252,7 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
// approximation of the derivatives.
const int n = 6;
const int np = phase_usage_.num_phases;
std::vector<int> pvtRegionIdx(n, 0);
// the tolerance for acceptable difference in values
const double reltol = 1e-9;
@@ -306,13 +297,13 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
}
testmu(reltol, n, np, p, r,z, props_, condition);
testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition);
testb(reltol,n,np,p,r,z,props_,condition);
testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition);
testrsSat(reltol,n,np,p,props_);
testrsSat(reltol,n,np,pvtRegionIdx,p,props_);
testrvSat(reltol,n,np,p,props_);
testrvSat(reltol,n,np,pvtRegionIdx,p,props_);
}
@@ -337,6 +328,7 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
// approximation of the derivatives.
const int n = 6;
const int np = phase_usage_.num_phases;
std::vector<int> pvtRegionIdx(n, 0);
// the tolerance for acceptable difference in values
const double reltol = 1e-9;
@@ -381,12 +373,12 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
}
testmu(reltol, n, np, p, r,z, props_, condition);
testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition);
testb(reltol,n,np,p,r,z,props_,condition);
testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition);
testrsSat(reltol,n,np,p,props_);
testrsSat(reltol,n,np,pvtRegionIdx,p,props_);
testrvSat(reltol,n,np,p,props_);
testrvSat(reltol,n,np,pvtRegionIdx,p,props_);
}