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

Merge pull request #688 from andlaus/temperature_dependent_PVT

Browse Source 
PVT properties: allow them to be temperature dependent
This commit is contained in:
Atgeirr Flø Rasmussen 2014-12-02 08:48:57 +01:00
commit fdad6860cb
35 changed files with 370 additions and 180 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
opm/core/pressure/CompressibleTpfa.cpp
View File

@ -237,7 +237,7 @@ namespace Opm
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w + 1]; ++j) { for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w + 1]; ++j) {
const int cell = wells_->well_cells[j]; const int cell = wells_->well_cells[j];
const double cell_depth = grid_.cell_centroids[dim * cell + dim - 1]; 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_.matrix(1, &state.pressure()[cell], &state.temperature()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0);
props_.density(1, &A[0], &cell, &rho[0]); props_.density(1, &A[0], &cell, &rho[0]);
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
const double s_phase = state.saturation()[np*cell + phase]; const double s_phase = state.saturation()[np*cell + phase];
@ -309,13 +309,14 @@ namespace Opm
const int nc = grid_.number_of_cells; const int nc = grid_.number_of_cells;
const int np = props_.numPhases(); const int np = props_.numPhases();
const double* cell_p = &state.pressure()[0]; const double* cell_p = &state.pressure()[0];
const double* cell_T = &state.temperature()[0];
const double* cell_z = &state.surfacevol()[0]; const double* cell_z = &state.surfacevol()[0];
const double* cell_s = &state.saturation()[0]; const double* cell_s = &state.saturation()[0];
cell_A_.resize(nc*np*np); cell_A_.resize(nc*np*np);
cell_dA_.resize(nc*np*np); cell_dA_.resize(nc*np*np);
props_.matrix(nc, cell_p, cell_z, &allcells_[0], &cell_A_[0], &cell_dA_[0]); props_.matrix(nc, cell_p, cell_T, cell_z, &allcells_[0], &cell_A_[0], &cell_dA_[0]);
cell_viscosity_.resize(nc*np); cell_viscosity_.resize(nc*np);
props_.viscosity(nc, cell_p, cell_z, &allcells_[0], &cell_viscosity_[0], 0); props_.viscosity(nc, cell_p, cell_T, cell_z, &allcells_[0], &cell_viscosity_[0], 0);
cell_phasemob_.resize(nc*np); cell_phasemob_.resize(nc*np);
props_.relperm(nc, cell_s, &allcells_[0], &cell_phasemob_[0], 0); props_.relperm(nc, cell_s, &allcells_[0], &cell_phasemob_[0], 0);
std::transform(cell_phasemob_.begin(), cell_phasemob_.end(), std::transform(cell_phasemob_.begin(), cell_phasemob_.end(),
@ -481,13 +482,14 @@ namespace Opm
} else { } else {
const double bhp = well_state.bhp()[w]; const double bhp = well_state.bhp()[w];
double perf_p = bhp + wellperf_wdp_[j]; double perf_p = bhp + wellperf_wdp_[j];
const double perf_T = well_state.temperature()[w];
// Hack warning: comp_frac is used as a component // Hack warning: comp_frac is used as a component
// surface-volume variable in calls to matrix() and // surface-volume variable in calls to matrix() and
// viscosity(), but as a saturation in the call to // viscosity(), but as a saturation in the call to
// relperm(). This is probably ok as long as injectors // relperm(). This is probably ok as long as injectors
// only inject pure fluids. // only inject pure fluids.
props_.matrix(1, &perf_p, comp_frac, &c, wpA, NULL); props_.matrix(1, &perf_p, &perf_T, comp_frac, &c, wpA, NULL);
props_.viscosity(1, &perf_p, comp_frac, &c, &mu[0], NULL); props_.viscosity(1, &perf_p, &perf_T, comp_frac, &c, &mu[0], NULL);
assert(std::fabs(std::accumulate(comp_frac, comp_frac + np, 0.0) - 1.0) < 1e-6); assert(std::fabs(std::accumulate(comp_frac, comp_frac + np, 0.0) - 1.0) < 1e-6);
props_.relperm (1, comp_frac, &c, wpM , NULL); props_.relperm (1, comp_frac, &c, wpM , NULL);
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {

10
opm/core/props/BlackoilPropertiesBasic.cpp
View File

@ -91,6 +91,7 @@ namespace Opm
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] p Array of n pressure values. /// \param[in] p Array of n pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] z Array of nP surface volume values. /// \param[in] z Array of nP surface volume values.
/// \param[in] cells Array of n cell indices to be associated with the p and z values. /// \param[in] cells Array of n cell indices to be associated with the p and z values.
/// \param[out] mu Array of nP viscosity values, array must be valid before calling. /// \param[out] mu Array of nP viscosity values, array must be valid before calling.
@ -98,6 +99,7 @@ namespace Opm
/// array must be valid before calling. /// array must be valid before calling.
void BlackoilPropertiesBasic::viscosity(const int n, void BlackoilPropertiesBasic::viscosity(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
const int* /*cells*/, const int* /*cells*/,
double* mu, double* mu,
@ -106,12 +108,13 @@ namespace Opm
if (dmudp) { if (dmudp) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesBasic::viscosity() -- derivatives of viscosity not yet implemented."); OPM_THROW(std::runtime_error, "BlackoilPropertiesBasic::viscosity() -- derivatives of viscosity not yet implemented.");
} else { } else {
pvt_.mu(n, p, z, mu); pvt_.mu(n, p, T, z, mu);
} }
} }
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] p Array of n pressure values. /// \param[in] p Array of n pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] z Array of nP surface volume values. /// \param[in] z Array of nP surface volume values.
/// \param[in] cells Array of n cell indices to be associated with the p and z values. /// \param[in] cells Array of n cell indices to be associated with the p and z values.
/// \param[out] A Array of nP^2 values, array must be valid before calling. /// \param[out] A Array of nP^2 values, array must be valid before calling.
@ -121,7 +124,8 @@ namespace Opm
/// array must be valid before calling. The matrices are output /// array must be valid before calling. The matrices are output
/// in Fortran order. /// in Fortran order.
void BlackoilPropertiesBasic::matrix(const int n, void BlackoilPropertiesBasic::matrix(const int n,
const double* /*p*/, const double* p,
const double* T,
const double* /*z*/, const double* /*z*/,
const int* /*cells*/, const int* /*cells*/,
double* A, double* A,
@ -130,7 +134,7 @@ namespace Opm
const int np = numPhases(); const int np = numPhases();
assert(np <= 2); assert(np <= 2);
double B[2]; // Must be enough since component classes do not handle more than 2. double B[2]; // Must be enough since component classes do not handle more than 2.
pvt_.B(1, 0, 0, B); pvt_.B(1, p, T, 0, B);
// Compute A matrix // Compute A matrix
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {

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

@ -83,6 +83,7 @@ namespace Opm
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] p Array of n pressure values. /// \param[in] p Array of n pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] z Array of nP surface volume values. /// \param[in] z Array of nP surface volume values.
/// \param[in] cells Array of n cell indices to be associated with the p and z values. /// \param[in] cells Array of n cell indices to be associated with the p and z values.
/// \param[out] mu Array of nP viscosity values, array must be valid before calling. /// \param[out] mu Array of nP viscosity values, array must be valid before calling.
@ -90,6 +91,7 @@ namespace Opm
/// array must be valid before calling. /// array must be valid before calling.
virtual void viscosity(const int n, virtual void viscosity(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
const int* cells, const int* cells,
double* mu, double* mu,
@ -97,6 +99,7 @@ namespace Opm
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] p Array of n pressure values. /// \param[in] p Array of n pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] z Array of nP surface volume values. /// \param[in] z Array of nP surface volume values.
/// \param[in] cells Array of n cell indices to be associated with the p and z values. /// \param[in] cells Array of n cell indices to be associated with the p and z values.
/// \param[out] A Array of nP^2 values, array must be valid before calling. /// \param[out] A Array of nP^2 values, array must be valid before calling.
@ -107,6 +110,7 @@ namespace Opm
/// in Fortran order. /// in Fortran order.
virtual void matrix(const int n, virtual void matrix(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
const int* cells, const int* cells,
double* A, double* A,

10
opm/core/props/BlackoilPropertiesFromDeck.cpp
View File

@ -90,6 +90,7 @@ namespace Opm
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] p Array of n pressure values. /// \param[in] p Array of n pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] z Array of nP surface volume values. /// \param[in] z Array of nP surface volume values.
/// \param[in] cells Array of n cell indices to be associated with the p and z values. /// \param[in] cells Array of n cell indices to be associated with the p and z values.
/// \param[out] mu Array of nP viscosity values, array must be valid before calling. /// \param[out] mu Array of nP viscosity values, array must be valid before calling.
@ -97,6 +98,7 @@ namespace Opm
/// array must be valid before calling. /// array must be valid before calling.
void BlackoilPropertiesFromDeck::viscosity(const int n, void BlackoilPropertiesFromDeck::viscosity(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
const int* cells, const int* cells,
double* mu, double* mu,
@ -111,12 +113,13 @@ namespace Opm
for (int i = 0; i < n; ++ i) for (int i = 0; i < n; ++ i)
pvtTableIdx[i] = cellPvtTableIdx[cells[i]]; pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
pvt_.mu(n, &pvtTableIdx[0], p, z, mu); pvt_.mu(n, &pvtTableIdx[0], p, T, z, mu);
} }
} }
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] p Array of n pressure values. /// \param[in] p Array of n pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] z Array of nP surface volume values. /// \param[in] z Array of nP surface volume values.
/// \param[in] cells Array of n cell indices to be associated with the p and z values. /// \param[in] cells Array of n cell indices to be associated with the p and z values.
/// \param[out] A Array of nP^2 values, array must be valid before calling. /// \param[out] A Array of nP^2 values, array must be valid before calling.
@ -127,6 +130,7 @@ namespace Opm
/// in Fortran order. /// in Fortran order.
void BlackoilPropertiesFromDeck::matrix(const int n, void BlackoilPropertiesFromDeck::matrix(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
const int* cells, const int* cells,
double* A, double* A,
@ -144,10 +148,10 @@ namespace Opm
if (dAdp) { if (dAdp) {
dB_.resize(n*np); dB_.resize(n*np);
dR_.resize(n*np); dR_.resize(n*np);
pvt_.dBdp(n, &pvtTableIdx[0], p, z, &B_[0], &dB_[0]); pvt_.dBdp(n, &pvtTableIdx[0], p, T, z, &B_[0], &dB_[0]);
pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]); pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]);
} else { } else {
pvt_.B(n, &pvtTableIdx[0], p, z, &B_[0]); pvt_.B(n, &pvtTableIdx[0], p, T, z, &B_[0]);
pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]); pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]);
} }
const int* phase_pos = pvt_.phasePosition(); const int* phase_pos = pvt_.phasePosition();

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

@ -125,6 +125,7 @@ namespace Opm
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] p Array of n pressure values. /// \param[in] p Array of n pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] z Array of nP surface volume values. /// \param[in] z Array of nP surface volume values.
/// \param[in] cells Array of n cell indices to be associated with the p and z values. /// \param[in] cells Array of n cell indices to be associated with the p and z values.
/// \param[out] mu Array of nP viscosity values, array must be valid before calling. /// \param[out] mu Array of nP viscosity values, array must be valid before calling.
@ -132,6 +133,7 @@ namespace Opm
/// array must be valid before calling. /// array must be valid before calling.
virtual void viscosity(const int n, virtual void viscosity(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
const int* cells, const int* cells,
double* mu, double* mu,
@ -139,6 +141,7 @@ namespace Opm
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] p Array of n pressure values. /// \param[in] p Array of n pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] z Array of nP surface volume values. /// \param[in] z Array of nP surface volume values.
/// \param[in] cells Array of n cell indices to be associated with the p and z values. /// \param[in] cells Array of n cell indices to be associated with the p and z values.
/// \param[out] A Array of nP^2 values, array must be valid before calling. /// \param[out] A Array of nP^2 values, array must be valid before calling.
@ -149,6 +152,7 @@ namespace Opm
/// in Fortran order. /// in Fortran order.
virtual void matrix(const int n, virtual void matrix(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
const int* cells, const int* cells,
double* A, double* A,

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

@ -70,6 +70,7 @@ namespace Opm
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] p Array of n pressure values. /// \param[in] p Array of n pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] z Array of nP surface volume values. /// \param[in] z Array of nP surface volume values.
/// \param[in] cells Array of n cell indices to be associated with the p and z values. /// \param[in] cells Array of n cell indices to be associated with the p and z values.
/// \param[out] mu Array of nP viscosity values, array must be valid before calling. /// \param[out] mu Array of nP viscosity values, array must be valid before calling.
@ -77,6 +78,7 @@ namespace Opm
/// array must be valid before calling. /// array must be valid before calling.
virtual void viscosity(const int n, virtual void viscosity(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
const int* cells, const int* cells,
double* mu, double* mu,
@ -84,6 +86,7 @@ namespace Opm
/// \param[in] n Number of data points. /// \param[in] n Number of data points.
/// \param[in] p Array of n pressure values. /// \param[in] p Array of n pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] z Array of nP surface volume values. /// \param[in] z Array of nP surface volume values.
/// \param[in] cells Array of n cell indices to be associated with the p and z values. /// \param[in] cells Array of n cell indices to be associated with the p and z values.
/// \param[out] A Array of nP^2 values, array must be valid before calling. /// \param[out] A Array of nP^2 values, array must be valid before calling.
@ -94,6 +97,7 @@ namespace Opm
/// in Fortran order. /// in Fortran order.
virtual void matrix(const int n, virtual void matrix(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
const int* cells, const int* cells,
double* A, double* A,

4
opm/core/props/IncompPropertiesBasic.cpp
View File

@ -44,7 +44,7 @@ namespace Opm
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
} }
viscosity_.resize(pvt_.numPhases()); viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]); pvt_.mu(1, 0, 0, 0, &viscosity_[0]);
} }
IncompPropertiesBasic::IncompPropertiesBasic(const int num_phases, IncompPropertiesBasic::IncompPropertiesBasic(const int num_phases,
@ -64,7 +64,7 @@ namespace Opm
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
} }
viscosity_.resize(pvt_.numPhases()); viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]); pvt_.mu(1, 0, 0, 0, &viscosity_[0]);
} }
IncompPropertiesBasic::~IncompPropertiesBasic() IncompPropertiesBasic::~IncompPropertiesBasic()

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

@ -161,12 +161,13 @@ namespace Opm
void BlackoilPvtProperties::mu(const int n, void BlackoilPvtProperties::mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_mu) const double* output_mu) const
{ {
data1_.resize(n); data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->mu(n, pvtTableIdx, p, z, &data1_[0]); props_[phase]->mu(n, pvtTableIdx, p, T, z, &data1_[0]);
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
output_mu[phase_usage_.num_phases*i + phase] = data1_[i]; output_mu[phase_usage_.num_phases*i + phase] = data1_[i];
@ -177,12 +178,13 @@ namespace Opm
void BlackoilPvtProperties::B(const int n, void BlackoilPvtProperties::B(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B) const double* output_B) const
{ {
data1_.resize(n); data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->B(n, pvtTableIdx, p, z, &data1_[0]); props_[phase]->B(n, pvtTableIdx, p, T, z, &data1_[0]);
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i]; output_B[phase_usage_.num_phases*i + phase] = data1_[i];
@ -193,6 +195,7 @@ namespace Opm
void BlackoilPvtProperties::dBdp(const int n, void BlackoilPvtProperties::dBdp(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B, double* output_B,
double* output_dBdp) const double* output_dBdp) const
@ -200,7 +203,7 @@ namespace Opm
data1_.resize(n); data1_.resize(n);
data2_.resize(n); data2_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->dBdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]); props_[phase]->dBdp(n, pvtTableIdx, p, T, z, &data1_[0], &data2_[0]);
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i]; output_B[phase_usage_.num_phases*i + phase] = data1_[i];

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

@ -77,24 +77,27 @@ namespace Opm
/// \return Array of size numPhases(). /// \return Array of size numPhases().
const double* surfaceDensities(int regionIdx = 0) const; const double* surfaceDensities(int regionIdx = 0) const;
/// Viscosity as a function of p and z. /// Viscosity as a function of p, T and z.
void mu(const int n, void mu(const int n,
const int *pvtTableIdx, const int *pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_mu) const; double* output_mu) const;
/// Formation volume factor as a function of p and z. /// Formation volume factor as a function of p, T and z.
void B(const int n, void B(const int n,
const int *pvtTableIdx, const int *pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B) const; double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z. /// Formation volume factor and p-derivative as functions of p, T and z.
void dBdp(const int n, void dBdp(const int n,
const int *pvtTableIdx, const int *pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B, double* output_B,
double* output_dBdp) const; double* output_dBdp) const;

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

@ -109,6 +109,7 @@ namespace Opm
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*z*/, const double* /*z*/,
double* output_mu) const double* output_mu) const
{ {
@ -124,6 +125,7 @@ namespace Opm
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
@ -144,6 +146,7 @@ namespace Opm
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
const PhasePresence* /*cond*/, const PhasePresence* /*cond*/,
double* output_mu, double* output_mu,
@ -165,6 +168,7 @@ namespace Opm
virtual void B(const int n, virtual void B(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*z*/, const double* /*z*/,
double* output_B) const double* output_B) const
{ {
@ -180,6 +184,7 @@ namespace Opm
virtual void dBdp(const int n, virtual void dBdp(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*z*/, const double* /*z*/,
double* output_B, double* output_B,
double* output_dBdp) const double* output_dBdp) const
@ -197,6 +202,7 @@ namespace Opm
virtual void b(const int n, virtual void b(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
double* output_b, double* output_b,
double* output_dbdp, double* output_dbdp,
@ -220,6 +226,7 @@ namespace Opm
virtual void b(const int n, virtual void b(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
const PhasePresence* /*cond*/, const PhasePresence* /*cond*/,
double* output_b, double* output_b,

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

@ -112,7 +112,8 @@ namespace Opm
void PvtDead::mu(const int n, void PvtDead::mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*z*/, const double* /*z*/,
double* output_mu) const double* output_mu) const
{ {
@ -127,7 +128,8 @@ namespace Opm
void PvtDead::mu(const int n, void PvtDead::mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
@ -148,7 +150,8 @@ namespace Opm
void PvtDead::mu(const int n, void PvtDead::mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
const PhasePresence* /*cond*/, const PhasePresence* /*cond*/,
double* output_mu, double* output_mu,
@ -171,7 +174,8 @@ namespace Opm
void PvtDead::B(const int n, void PvtDead::B(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*z*/, const double* /*z*/,
double* output_B) const double* output_B) const
{ {
@ -185,12 +189,13 @@ namespace Opm
void PvtDead::dBdp(const int n, void PvtDead::dBdp(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* /*z*/, const double* /*z*/,
double* output_B, double* output_B,
double* output_dBdp) const double* output_dBdp) const
{ {
B(n, pvtTableIdx, p, 0, output_B); B(n, pvtTableIdx, p, T, 0, output_B);
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i); int regionIdx = getTableIndex_(pvtTableIdx, i);
@ -201,7 +206,8 @@ namespace Opm
void PvtDead::b(const int n, void PvtDead::b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
double* output_b, double* output_b,
double* output_dbdp, double* output_dbdp,
@ -222,7 +228,8 @@ namespace Opm
void PvtDead::b(const int n, void PvtDead::b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
const PhasePresence* /*cond*/, const PhasePresence* /*cond*/,
double* output_b, double* output_b,

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

@ -50,64 +50,71 @@ namespace Opm
void initFromGas(const std::vector<Opm::PvdgTable>& pvdgTables); void initFromGas(const std::vector<Opm::PvdgTable>& pvdgTables);
virtual ~PvtDead(); virtual ~PvtDead();
/// Viscosity as a function of p and z. /// Viscosity as a function of p, T and z.
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_mu) const; double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r. /// Viscosity and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p). /// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
double* output_dmudr) const; double* output_dmudr) const;
/// Viscosity as a function of p and r. /// Viscosity as a function of p, T and r.
/// State condition determined by 'cond'. /// State condition determined by 'cond'.
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
double* output_dmudr) const; double* output_dmudr) const;
/// Formation volume factor as a function of p and z. /// Formation volume factor as a function of p, T and z.
virtual void B(const int n, virtual void B(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B) const; double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z. /// Formation volume factor and p-derivative as functions of p, T and z.
virtual void dBdp(const int n, virtual void dBdp(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B, double* output_B,
double* output_dBdp) const; double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p). /// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n, virtual void b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
double* output_b, double* output_b,
double* output_dbdp, double* output_dbdp,
double* output_dbdr) const; double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// State condition determined by 'cond'. /// State condition determined by 'cond'.
virtual void b(const int n, virtual void b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_b, double* output_b,

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

@ -105,7 +105,8 @@ namespace Opm
void PvtDeadSpline::mu(const int n, void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*z*/, const double* /*z*/,
double* output_mu) const double* output_mu) const
{ {
@ -119,6 +120,7 @@ namespace Opm
void PvtDeadSpline::mu(const int n, void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
@ -136,6 +138,7 @@ namespace Opm
void PvtDeadSpline::mu(const int n, void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
const PhasePresence* /*cond*/, const PhasePresence* /*cond*/,
double* output_mu, double* output_mu,
@ -154,7 +157,8 @@ namespace Opm
void PvtDeadSpline::B(const int n, void PvtDeadSpline::B(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* /*z*/, const double* /*z*/,
double* output_B) const double* output_B) const
{ {
@ -168,11 +172,12 @@ namespace Opm
void PvtDeadSpline::dBdp(const int n, void PvtDeadSpline::dBdp(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* /*z*/, const double* /*z*/,
double* output_B, double* output_B,
double* output_dBdp) const double* output_dBdp) const
{ {
B(n, pvtTableIdx, p, 0, output_B); B(n, pvtTableIdx, p, T, 0, output_B);
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i); int regionIdx = getTableIndex_(pvtTableIdx, i);
@ -183,7 +188,8 @@ namespace Opm
void PvtDeadSpline::b(const int n, void PvtDeadSpline::b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* /*r*/, const double* /*r*/,
double* output_b, double* output_b,
double* output_dbdp, double* output_dbdp,
@ -201,6 +207,7 @@ namespace Opm
void PvtDeadSpline::b(const int n, void PvtDeadSpline::b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* /*r*/, const double* /*r*/,
const PhasePresence* /*cond*/, const PhasePresence* /*cond*/,
double* output_b, double* output_b,

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

@ -49,64 +49,71 @@ namespace Opm
virtual ~PvtDeadSpline(); virtual ~PvtDeadSpline();
/// Viscosity as a function of p and z. /// Viscosity as a function of p, T and z.
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_mu) const; double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r. /// Viscosity and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p). /// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
double* output_dmudr) const; double* output_dmudr) const;
/// Viscosity as a function of p and r. /// Viscosity as a function of p, T and r.
/// State condition determined by 'cond'. /// State condition determined by 'cond'.
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
double* output_dmudr) const; double* output_dmudr) const;
/// Formation volume factor as a function of p and z. /// Formation volume factor as a function of p, T and z.
virtual void B(const int n, virtual void B(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B) const; double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z. /// Formation volume factor and p-derivative as functions of p, T and z.
virtual void dBdp(const int n, virtual void dBdp(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B, double* output_B,
double* output_dBdp) const; double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p). /// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n, virtual void b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
double* output_b, double* output_b,
double* output_dbdp, double* output_dbdp,
double* output_dbdr) const; double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// State condition determined by 'cond'. /// State condition determined by 'cond'.
virtual void b(const int n, virtual void b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_b, double* output_b,

31
opm/core/props/pvt/PvtInterface.hpp
View File

@ -42,8 +42,8 @@ namespace Opm
/// arbitrary two-phase and three-phase situations. /// arbitrary two-phase and three-phase situations.
void setPhaseConfiguration(const int num_phases, const int* phase_pos); void setPhaseConfiguration(const int num_phases, const int* phase_pos);
/// The PVT properties can either be given as a function of pressure (p) and surface volume (z) /// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z)
/// or pressure (p) and gas resolution factor (r). /// or pressure (p), temperature (T) and gas resolution factor (r).
/// For all the virtual methods, the following apply: /// For all the virtual methods, the following apply:
/// - pvtRegionIdx is an array of size n and represents the /// - pvtRegionIdx is an array of size n and represents the
/// index of the PVT table which should be used to calculate /// index of the PVT table which should be used to calculate
@ -54,38 +54,42 @@ namespace Opm
/// - Output arrays shall be of size n, and must be valid before /// - Output arrays shall be of size n, and must be valid before
/// calling the method. /// calling the method.
/// Viscosity as a function of p and z. /// Viscosity as a function of p, T and z.
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_mu) const = 0; double* output_mu) const = 0;
/// Viscosity as a function of p and r. /// Viscosity as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p). /// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
double* output_dmudr) const = 0; double* output_dmudr) const = 0;
/// Viscosity as a function of p and r. /// Viscosity as a function of p, T and r.
/// State condition determined by 'cond'. /// State condition determined by 'cond'.
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
double* output_dmudr) const = 0; double* output_dmudr) const = 0;
/// Formation volume factor as a function of p and z. /// Formation volume factor as a function of p, T and z.
virtual void B(const int n, virtual void B(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B) const = 0; double* output_B) const = 0;
@ -93,25 +97,28 @@ namespace Opm
virtual void dBdp(const int n, virtual void dBdp(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B, double* output_B,
double* output_dBdp) const = 0; double* output_dBdp) const = 0;
/// The inverse of the volume factor b = 1 / B as a function of p and r. /// The inverse of the volume factor b = 1 / B as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p). /// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n, virtual void b(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
double* output_b, double* output_b,
double* output_dbdp, double* output_dbdp,
double* output_dpdr) const = 0; double* output_dpdr) const = 0;
/// The inverse of the volume factor b = 1 / B as a function of p and r. /// The inverse of the volume factor b = 1 / B as a function of p, T and r.
/// State condition determined by 'cond'. /// State condition determined by 'cond'.
virtual void b(const int n, virtual void b(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_b, double* output_b,

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

@ -104,6 +104,7 @@ namespace Opm
void PvtLiveGas::mu(const int n, void PvtLiveGas::mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* z, const double* z,
double* output_mu) const double* output_mu) const
{ {
@ -116,10 +117,11 @@ namespace Opm
} }
} }
/// Viscosity and its derivatives as a function of p and r. /// Viscosity and its p and r derivatives as a function of p, T and r.
void PvtLiveGas::mu(const int /*n*/, void PvtLiveGas::mu(const int /*n*/,
const int* /*pvtRegionIdx*/, const int* /*pvtRegionIdx*/,
const double* /*p*/, const double* /*p*/,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
double* /*output_mu*/, double* /*output_mu*/,
double* /*output_dmudp*/, double* /*output_dmudp*/,
@ -128,10 +130,11 @@ namespace Opm
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented"); OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
} }
/// Viscosity and its derivatives as a function of p and r. /// Viscosity and its p and r derivatives as a function of p, T and r.
void PvtLiveGas::mu(const int n, void PvtLiveGas::mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_mu, double* output_mu,
@ -164,10 +167,11 @@ namespace Opm
} }
/// Formation volume factor as a function of p and z. /// Formation volume factor as a function of p, T and z.
void PvtLiveGas::B(const int n, void PvtLiveGas::B(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* z, const double* z,
double* output_B) const double* output_B) const
{ {
@ -181,8 +185,9 @@ namespace Opm
/// Formation volume factor and p-derivative as functions of p and z. /// Formation volume factor and p-derivative as functions of p and z.
void PvtLiveGas::dBdp(const int n, void PvtLiveGas::dBdp(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* z, const double* z,
double* output_B, double* output_B,
double* output_dBdp) const double* output_dBdp) const
@ -193,10 +198,11 @@ namespace Opm
} }
} }
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
void PvtLiveGas::b(const int /*n*/, void PvtLiveGas::b(const int /*n*/,
const int* /*pvtRegionIdx*/, const int* /*pvtRegionIdx*/,
const double* /*p*/, const double* /*p*/,
const double* /*T*/,
const double* /*r*/, const double* /*r*/,
double* /*output_b*/, double* /*output_b*/,
double* /*output_dbdp*/, double* /*output_dbdp*/,
@ -206,10 +212,11 @@ namespace Opm
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented"); 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. /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
void PvtLiveGas::b(const int n, void PvtLiveGas::b(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* /*T*/,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_b, double* output_b,

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

@ -29,8 +29,8 @@
namespace Opm namespace Opm
{ {
/// Class for miscible wet gas (with vaporized oil in vapour phase). /// Class for miscible wet gas (with vaporized oil in vapour phase).
/// The PVT properties can either be given as a function of pressure (p) and surface volume (z) /// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z)
/// or pressure (p) and gas resolution factor (r). /// or pressure (p), temperature (T) and gas resolution factor (r).
/// For all the virtual methods, the following apply: p, r and z /// 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. /// 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 /// Output arrays shall be of size n, and must be valid before
@ -41,64 +41,71 @@ namespace Opm
PvtLiveGas(const std::vector<Opm::PvtgTable>& pvtgTables); PvtLiveGas(const std::vector<Opm::PvtgTable>& pvtgTables);
virtual ~PvtLiveGas(); virtual ~PvtLiveGas();
/// Viscosity as a function of p and z. /// Viscosity as a function of p, T and z.
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_mu) const; double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r. /// Viscosity and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p). /// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
double* output_dmudr) const; double* output_dmudr) const;
/// Viscosity as a function of p and r. /// Viscosity as a function of p, T and r.
/// State condition determined by 'cond'. /// State condition determined by 'cond'.
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
double* output_dmudr) const; double* output_dmudr) const;
/// Formation volume factor as a function of p and z. /// Formation volume factor as a function of p, T and z.
virtual void B(const int n, virtual void B(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B) const; double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z. /// Formation volume factor and p-derivative as functions of p, T and z.
virtual void dBdp(const int n, virtual void dBdp(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B, double* output_B,
double* output_dBdp) const; double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p). /// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n, virtual void b(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
double* output_b, double* output_b,
double* output_dbdp, double* output_dbdp,
double* output_dbdr) const; double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// State condition determined by 'cond'. /// State condition determined by 'cond'.
virtual void b(const int n, virtual void b(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_b, double* output_b,
@ -107,7 +114,7 @@ namespace Opm
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p. /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p, T.
virtual void rsSat(const int n, virtual void rsSat(const int n,
const int* pvtRegionIdx, const int* pvtRegionIdx,
const double* p, const double* p,

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

@ -132,10 +132,11 @@ namespace Opm
} }
/// Viscosity as a function of p and z. /// Viscosity as a function of p, T and z.
void PvtLiveOil::mu(const int n, void PvtLiveOil::mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* z, const double* z,
double* output_mu) const double* output_mu) const
{ {
@ -150,10 +151,11 @@ namespace Opm
} }
} }
/// Viscosity and its derivatives as a function of p and r. /// Viscosity and its p and r derivatives as a function of p, T and r.
void PvtLiveOil::mu(const int n, void PvtLiveOil::mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* r, const double* r,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
@ -184,10 +186,11 @@ namespace Opm
} }
} }
/// Viscosity and its derivatives as a function of p and r. /// Viscosity and its p and r derivatives as a function of p, T and r.
void PvtLiveOil::mu(const int n, void PvtLiveOil::mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_mu, double* output_mu,
@ -221,10 +224,11 @@ namespace Opm
} }
/// Formation volume factor as a function of p and z. /// Formation volume factor as a function of p, T and z.
void PvtLiveOil::B(const int n, void PvtLiveOil::B(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* z, const double* z,
double* output_B) const double* output_B) const
{ {
@ -238,10 +242,11 @@ namespace Opm
} }
/// Formation volume factor and p-derivative as functions of p and z. /// Formation volume factor and p-derivative as functions of p, T and z.
void PvtLiveOil::dBdp(const int n, void PvtLiveOil::dBdp(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* z, const double* z,
double* output_B, double* output_B,
double* output_dBdp) const double* output_dBdp) const
@ -255,8 +260,9 @@ namespace Opm
} }
void PvtLiveOil::b(const int n, void PvtLiveOil::b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* r, const double* r,
double* output_b, double* output_b,
double* output_dbdp, double* output_dbdp,
@ -275,8 +281,9 @@ namespace Opm
} }
void PvtLiveOil::b(const int n, void PvtLiveOil::b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* /*T*/,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_b, double* output_b,

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

@ -30,8 +30,8 @@
namespace Opm namespace Opm
{ {
/// Class for miscible live oil (with dissolved gas in liquid phase). /// Class for miscible live oil (with dissolved gas in liquid phase).
/// The PVT properties can either be given as a function of pressure (p) and surface volume (z) /// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z)
/// or pressure (p) and gas resolution factor (r). /// or pressure (p), temperature (T) and gas resolution factor (r).
/// For all the virtual methods, the following apply: p, r and z /// 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. /// 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 /// Output arrays shall be of size n, and must be valid before
@ -42,64 +42,71 @@ namespace Opm
PvtLiveOil(const std::vector<Opm::PvtoTable>& pvtoTables); PvtLiveOil(const std::vector<Opm::PvtoTable>& pvtoTables);
virtual ~PvtLiveOil(); virtual ~PvtLiveOil();
/// Viscosity as a function of p and z. /// Viscosity as a function of p, T and z.
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_mu) const; double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r. /// Viscosity and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p). /// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
double* output_dmudr) const; double* output_dmudr) const;
/// Viscosity as a function of p and r. /// Viscosity as a function of p, T and r.
/// State condition determined by 'cond'. /// State condition determined by 'cond'.
virtual void mu(const int n, virtual void mu(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_mu, double* output_mu,
double* output_dmudp, double* output_dmudp,
double* output_dmudr) const; double* output_dmudr) const;
/// Formation volume factor as a function of p and z. /// Formation volume factor as a function of p, T and z.
virtual void B(const int n, virtual void B(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B) const; double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z. /// Formation volume factor and p-derivative as functions of p, T and z.
virtual void dBdp(const int n, virtual void dBdp(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B, double* output_B,
double* output_dBdp) const; double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p). /// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n, virtual void b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
double* output_b, double* output_b,
double* output_dbdp, double* output_dbdp,
double* output_dbdr) const; double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p, T and r.
/// State condition determined by 'cond'. /// State condition determined by 'cond'.
virtual void b(const int n, virtual void b(const int n,
const int* pvtTableIdx, const int* pvtTableIdx,
const double* p, const double* p,
const double* T,
const double* r, const double* r,
const PhasePresence* cond, const PhasePresence* cond,
double* output_b, double* output_b,

3
opm/core/props/pvt/PvtPropertiesBasic.cpp
View File

@ -113,6 +113,7 @@ namespace Opm
void PvtPropertiesBasic::mu(const int n, void PvtPropertiesBasic::mu(const int n,
const double* /*p*/, const double* /*p*/,
const double* /*T*/,
const double* /*z*/, const double* /*z*/,
double* output_mu) const double* output_mu) const
{ {
@ -127,6 +128,7 @@ namespace Opm
void PvtPropertiesBasic::B(const int n, void PvtPropertiesBasic::B(const int n,
const double* /*p*/, const double* /*p*/,
const double* /*T*/,
const double* /*z*/, const double* /*z*/,
double* output_B) const double* output_B) const
{ {
@ -141,6 +143,7 @@ namespace Opm
void PvtPropertiesBasic::dBdp(const int n, void PvtPropertiesBasic::dBdp(const int n,
const double* /*p*/, const double* /*p*/,
const double* /*T*/,
const double* /*z*/, const double* /*z*/,
double* output_B, double* output_B,
double* output_dBdp) const double* output_dBdp) const

11
opm/core/props/pvt/PvtPropertiesBasic.hpp
View File

@ -29,7 +29,7 @@ namespace Opm
/// Class collecting simple pvt properties for 1-3 phases. /// Class collecting simple pvt properties for 1-3 phases.
/// All phases are incompressible and have constant viscosities. /// All phases are incompressible and have constant viscosities.
/// For all the methods, the following apply: p and z are unused. /// For all the methods, the following apply: p, T and z are unused.
/// Output arrays shall be of size n*numPhases(), and must be valid /// Output arrays shall be of size n*numPhases(), and must be valid
/// before calling the method. /// before calling the method.
/// NOTE: This class is intentionally similar to BlackoilPvtProperties. /// NOTE: This class is intentionally similar to BlackoilPvtProperties.
@ -62,21 +62,24 @@ namespace Opm
/// \return Array of size numPhases(). /// \return Array of size numPhases().
const double* surfaceDensities() const; const double* surfaceDensities() const;
/// Viscosity as a function of p and z. /// Viscosity as a function of p, T and z.
void mu(const int n, void mu(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_mu) const; double* output_mu) const;
/// Formation volume factor as a function of p and z. /// Formation volume factor as a function of p, T and z.
void B(const int n, void B(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B) const; double* output_B) const;
/// Formation volume factor and p-derivative as functions of p and z. /// Formation volume factor and p-derivative as functions of p, T and z.
void dBdp(const int n, void dBdp(const int n,
const double* p, const double* p,
const double* T,
const double* z, const double* z,
double* output_B, double* output_B,
double* output_dBdp) const; double* output_dBdp) const;

71
opm/core/simulator/EquilibrationHelpers.hpp
View File

@ -122,12 +122,15 @@ namespace Opm
* *
* \param[in] p Fluid pressure. * \param[in] p Fluid pressure.
* *
* \param[in] T Temperature.
*
* \param[in] z Surface volumes of all phases. * \param[in] z Surface volumes of all phases.
* *
* \return Phase densities at phase point. * \return Phase densities at phase point.
*/ */
std::vector<double> std::vector<double>
operator()(const double p, operator()(const double p,
const double T,
const std::vector<double>& z) const const std::vector<double>& z) const
{ {
const int np = props_.numPhases(); const int np = props_.numPhases();
@ -136,7 +139,7 @@ namespace Opm
assert (z.size() == std::vector<double>::size_type(np)); assert (z.size() == std::vector<double>::size_type(np));
double* dAdp = 0; double* dAdp = 0;
props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp); props_.matrix(1, &p, &T, &z[0], &c_[0], &A[0], dAdp);
std::vector<double> rho(np, 0.0); std::vector<double> rho(np, 0.0);
props_.density(1, &A[0], &c_[0], &rho[0]); props_.density(1, &A[0], &c_[0], &rho[0]);
@ -171,11 +174,15 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS * \param[in] press Pressure at which to calculate RS
* value. * value.
* *
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c * \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
virtual double operator()(const double depth, virtual double operator()(const double depth,
const double press, const double press,
const double temp,
const double sat = 0.0) const = 0; const double sat = 0.0) const = 0;
}; };
@ -194,6 +201,9 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS * \param[in] press Pressure at which to calculate RS
* value. * value.
* *
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c * \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. In "no mixing * depth and pressure @c press. In "no mixing
* policy", this is identically zero. * policy", this is identically zero.
@ -201,6 +211,7 @@ namespace Opm
double double
operator()(const double /* depth */, operator()(const double /* depth */,
const double /* press */, const double /* press */,
const double /* temp */,
const double /* sat */ = 0.0) const const double /* sat */ = 0.0) const
{ {
return 0.0; return 0.0;
@ -247,18 +258,22 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS * \param[in] press Pressure at which to calculate RS
* value. * value.
* *
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c * \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double double
operator()(const double depth, operator()(const double depth,
const double press, const double press,
const double temp,
const double sat_gas = 0.0) const const double sat_gas = 0.0) const
{ {
if (sat_gas > 0.0) { if (sat_gas > 0.0) {
return satRs(press); return satRs(press, temp);
} else { } else {
return std::min(satRs(press), linearInterpolation(depth_, rs_, depth)); return std::min(satRs(press, temp), linearInterpolation(depth_, rs_, depth));
} }
} }
@ -270,9 +285,9 @@ namespace Opm
double z_[BlackoilPhases::MaxNumPhases]; double z_[BlackoilPhases::MaxNumPhases];
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases]; mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRs(const double press) const double satRs(const double press, const double temp) const
{ {
props_.matrix(1, &press, z_, &cell_, A_, 0); props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rs/Bo is in the gas row and oil column of A_. // Rs/Bo is in the gas row and oil column of A_.
// 1/Bo is in the oil row and column. // 1/Bo is in the oil row and column.
// Recall also that it is stored in column-major order. // Recall also that it is stored in column-major order.
@ -323,18 +338,22 @@ namespace Opm
* \param[in] press Pressure at which to calculate RV * \param[in] press Pressure at which to calculate RV
* value. * value.
* *
* \param[in] temp Temperature at which to calculate RV
* value.
*
* \return Vaporized oil-gas ratio (RV) at depth @c * \return Vaporized oil-gas ratio (RV) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double double
operator()(const double depth, operator()(const double depth,
const double press, const double press,
const double temp,
const double sat_oil = 0.0 ) const const double sat_oil = 0.0 ) const
{ {
if (sat_oil > 0.0) { if (sat_oil > 0.0) {
return satRv(press); return satRv(press, temp);
} else { } else {
return std::min(satRv(press), linearInterpolation(depth_, rv_, depth)); return std::min(satRv(press, temp), linearInterpolation(depth_, rv_, depth));
} }
} }
@ -346,9 +365,9 @@ namespace Opm
double z_[BlackoilPhases::MaxNumPhases]; double z_[BlackoilPhases::MaxNumPhases];
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases]; mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRv(const double press) const double satRv(const double press, const double temp) const
{ {
props_.matrix(1, &press, z_, &cell_, A_, 0); props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rv/Bg is in the oil row and gas column of A_. // Rv/Bg is in the oil row and gas column of A_.
// 1/Bg is in the gas row and column. // 1/Bg is in the gas row and column.
// Recall also that it is stored in column-major order. // Recall also that it is stored in column-major order.
@ -382,15 +401,16 @@ namespace Opm
* \param[in] props property object * \param[in] props property object
* \param[in] cell any cell in the pvt region * \param[in] cell any cell in the pvt region
* \param[in] p_contact oil pressure at the contact * \param[in] p_contact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/ */
RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact) RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
: props_(props), cell_(cell) : props_(props), cell_(cell)
{ {
auto pu = props_.phaseUsage(); auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0); std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100; z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
rs_sat_contact_ = satRs(p_contact); rs_sat_contact_ = satRs(p_contact, T_contact);
} }
/** /**
@ -402,18 +422,22 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS * \param[in] press Pressure at which to calculate RS
* value. * value.
* *
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c * \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double double
operator()(const double /* depth */, operator()(const double /* depth */,
const double press, const double press,
const double temp,
const double sat_gas = 0.0) const const double sat_gas = 0.0) const
{ {
if (sat_gas > 0.0) { if (sat_gas > 0.0) {
return satRs(press); return satRs(press, temp);
} else { } else {
return std::min(satRs(press), rs_sat_contact_); return std::min(satRs(press, temp), rs_sat_contact_);
} }
} }
@ -424,9 +448,9 @@ namespace Opm
double rs_sat_contact_; double rs_sat_contact_;
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases]; mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRs(const double press) const double satRs(const double press, const double temp) const
{ {
props_.matrix(1, &press, z_, &cell_, A_, 0); props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rs/Bo is in the gas row and oil column of A_. // Rs/Bo is in the gas row and oil column of A_.
// 1/Bo is in the oil row and column. // 1/Bo is in the oil row and column.
// Recall also that it is stored in column-major order. // Recall also that it is stored in column-major order.
@ -460,15 +484,16 @@ namespace Opm
* \param[in] props property object * \param[in] props property object
* \param[in] cell any cell in the pvt region * \param[in] cell any cell in the pvt region
* \param[in] p_contact oil pressure at the contact * \param[in] p_contact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/ */
RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact) RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
: props_(props), cell_(cell) : props_(props), cell_(cell)
{ {
auto pu = props_.phaseUsage(); auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0); std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100; z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
rv_sat_contact_ = satRv(p_contact); rv_sat_contact_ = satRv(p_contact, T_contact);
} }
/** /**
@ -480,18 +505,22 @@ namespace Opm
* \param[in] press Pressure at which to calculate RV * \param[in] press Pressure at which to calculate RV
* value. * value.
* *
* \param[in] temp Temperature at which to calculate RV
* value.
*
* \return Dissolved oil-gas ratio (RV) at depth @c * \return Dissolved oil-gas ratio (RV) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double double
operator()(const double /*depth*/, operator()(const double /*depth*/,
const double press, const double press,
const double temp,
const double sat_oil = 0.0) const const double sat_oil = 0.0) const
{ {
if (sat_oil > 0.0) { if (sat_oil > 0.0) {
return satRv(press); return satRv(press, temp);
} else { } else {
return std::min(satRv(press), rv_sat_contact_); return std::min(satRv(press, temp), rv_sat_contact_);
} }
} }
@ -502,9 +531,9 @@ namespace Opm
double rv_sat_contact_; double rv_sat_contact_;
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases]; mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRv(const double press) const double satRv(const double press, const double temp) const
{ {
props_.matrix(1, &press, z_, &cell_, A_, 0); props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rv/Bg is in the oil row and gas column of A_. // Rv/Bg is in the oil row and gas column of A_.
// 1/Bg is in the gas row and column. // 1/Bg is in the gas row and column.
// Recall also that it is stored in column-major order. // Recall also that it is stored in column-major order.

4
opm/core/simulator/SimulatorCompressibleTwophase.cpp
View File

@ -358,7 +358,7 @@ namespace Opm
// Solve pressure equation. // Solve pressure equation.
if (check_well_controls_) { if (check_well_controls_) {
computeFractionalFlow(props_, allcells_, computeFractionalFlow(props_, allcells_,
state.pressure(), state.surfacevol(), state.saturation(), state.pressure(), state.temperature(), state.surfacevol(), state.saturation(),
fractional_flows); fractional_flows);
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
} }
@ -445,7 +445,7 @@ namespace Opm
double injected[2] = { 0.0 }; double injected[2] = { 0.0 };
double produced[2] = { 0.0 }; double produced[2] = { 0.0 };
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], &state.pressure()[0], tsolver_.solve(&state.faceflux()[0], &state.pressure()[0], &state.temperature()[0],
&initial_porevol[0], &porevol[0], &transport_src[0], stepsize, &initial_porevol[0], &porevol[0], &transport_src[0], stepsize,
state.saturation(), state.surfacevol()); state.saturation(), state.surfacevol());
double substep_injected[2] = { 0.0 }; double substep_injected[2] = { 0.0 };

2
opm/core/simulator/SimulatorState.cpp
View File

@ -12,6 +12,7 @@ SimulatorState::equals (const SimulatorState& other,
// if we use &=, then all the tests will be run regardless // if we use &=, then all the tests will be run regardless
equal = equal && vectorApproxEqual( pressure() , other.pressure() , epsilon); equal = equal && vectorApproxEqual( pressure() , other.pressure() , epsilon);
equal = equal && vectorApproxEqual( temperature() , other.temperature() , epsilon);
equal = equal && vectorApproxEqual( facepressure() , other.facepressure() , epsilon); equal = equal && vectorApproxEqual( facepressure() , other.facepressure() , epsilon);
equal = equal && vectorApproxEqual( faceflux() , other.faceflux() , epsilon); equal = equal && vectorApproxEqual( faceflux() , other.faceflux() , epsilon);
equal = equal && vectorApproxEqual( saturation() , other.saturation() , epsilon); equal = equal && vectorApproxEqual( saturation() , other.saturation() , epsilon);
@ -48,6 +49,7 @@ SimulatorState::init(int number_of_cells, int number_of_faces, int num_phases)
{ {
num_phases_ = num_phases; num_phases_ = num_phases;
press_.resize(number_of_cells, 0.0); press_.resize(number_of_cells, 0.0);
temp_.resize(number_of_cells, 273.15 + 20);
fpress_.resize(number_of_faces, 0.0); fpress_.resize(number_of_faces, 0.0);
flux_.resize(number_of_faces, 0.0); flux_.resize(number_of_faces, 0.0);
sat_.resize(num_phases_ * number_of_cells, 0.0); sat_.resize(num_phases_ * number_of_cells, 0.0);

3
opm/core/simulator/SimulatorState.hpp
View File

@ -38,11 +38,13 @@ namespace Opm
int numPhases() const { return num_phases_; } int numPhases() const { return num_phases_; }
std::vector<double>& pressure () { return press_ ; } std::vector<double>& pressure () { return press_ ; }
std::vector<double>& temperature () { return temp_ ; }
std::vector<double>& facepressure() { return fpress_; } std::vector<double>& facepressure() { return fpress_; }
std::vector<double>& faceflux () { return flux_ ; } std::vector<double>& faceflux () { return flux_ ; }
std::vector<double>& saturation () { return sat_ ; } std::vector<double>& saturation () { return sat_ ; }
const std::vector<double>& pressure () const { return press_ ; } const std::vector<double>& pressure () const { return press_ ; }
const std::vector<double>& temperature () const { return temp_ ; }
const std::vector<double>& facepressure() const { return fpress_; } const std::vector<double>& facepressure() const { return fpress_; }
const std::vector<double>& faceflux () const { return flux_ ; } const std::vector<double>& faceflux () const { return flux_ ; }
const std::vector<double>& saturation () const { return sat_ ; } const std::vector<double>& saturation () const { return sat_ ; }
@ -56,6 +58,7 @@ namespace Opm
private: private:
int num_phases_; int num_phases_;
std::vector<double> press_ ; std::vector<double> press_ ;
std::vector<double> temp_ ;
std::vector<double> fpress_; std::vector<double> fpress_;
std::vector<double> flux_ ; std::vector<double> flux_ ;
std::vector<double> sat_ ; std::vector<double> sat_ ;

6
opm/core/simulator/WellState.hpp
View File

@ -44,6 +44,7 @@ namespace Opm
const int nw = wells->number_of_wells; const int nw = wells->number_of_wells;
const int np = wells->number_of_phases; const int np = wells->number_of_phases;
bhp_.resize(nw); bhp_.resize(nw);
temperature_.resize(nw, 273.15 + 20); // standard temperature for now
wellrates_.resize(nw * np, 0.0); wellrates_.resize(nw * np, 0.0);
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER)); assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
@ -110,6 +111,10 @@ namespace Opm
std::vector<double>& bhp() { return bhp_; } std::vector<double>& bhp() { return bhp_; }
const std::vector<double>& bhp() const { return bhp_; } const std::vector<double>& bhp() const { return bhp_; }
/// One temperature per well.
std::vector<double>& temperature() { return temperature_; }
const std::vector<double>& temperature() const { return temperature_; }
/// One rate per well and phase. /// One rate per well and phase.
std::vector<double>& wellRates() { return wellrates_; } std::vector<double>& wellRates() { return wellrates_; }
const std::vector<double>& wellRates() const { return wellrates_; } const std::vector<double>& wellRates() const { return wellrates_; }
@ -124,6 +129,7 @@ namespace Opm
private: private:
std::vector<double> bhp_; std::vector<double> bhp_;
std::vector<double> temperature_;
std::vector<double> wellrates_; std::vector<double> wellrates_;
std::vector<double> perfrates_; std::vector<double> perfrates_;
std::vector<double> perfpress_; std::vector<double> perfpress_;

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

@ -171,6 +171,7 @@ namespace Opm
* \param[in] cells Range that spans the cells of the current * \param[in] cells Range that spans the cells of the current
* equilibration region. * equilibration region.
* \param[in] oil_pressure Oil pressure for each cell in range. * \param[in] oil_pressure Oil pressure for each cell in range.
* \param[in] temperature Temperature for each cell in range.
* \param[in] rs_func Rs as function of pressure and depth. * \param[in] rs_func Rs as function of pressure and depth.
* \return Rs values, one for each cell in the 'cells' range. * \return Rs values, one for each cell in the 'cells' range.
*/ */
@ -178,6 +179,7 @@ namespace Opm
std::vector<double> computeRs(const UnstructuredGrid& grid, std::vector<double> computeRs(const UnstructuredGrid& grid,
const CellRangeType& cells, const CellRangeType& cells,
const std::vector<double> oil_pressure, const std::vector<double> oil_pressure,
const std::vector<double>& temperature,
const Miscibility::RsFunction& rs_func, const Miscibility::RsFunction& rs_func,
const std::vector<double> gas_saturation); const std::vector<double> gas_saturation);
@ -298,7 +300,8 @@ namespace Opm
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold."); "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
} }
const double p_contact = rec[i].main.press; const double p_contact = rec[i].main.press;
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact)); const double T_contact = 273.15 + 20; // standard temperature for now
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact, T_contact));
} }
} }
} else { } else {
@ -329,7 +332,8 @@ namespace Opm
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold."); "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
} }
const double p_contact = rec[i].main.press + rec[i].goc.press; const double p_contact = rec[i].main.press + rec[i].goc.press;
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact)); const double T_contact = 273.15 + 20; // standard temperature for now
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact, T_contact));
} }
} }
} else { } else {
@ -399,6 +403,7 @@ namespace Opm
props.phaseUsage()); props.phaseUsage());
PVec press = phasePressures(G, eqreg, cells, grav); PVec press = phasePressures(G, eqreg, cells, grav);
const std::vector<double>& temp = temperature(G, eqreg, cells);
const PVec sat = phaseSaturations(G, eqreg, cells, props, swat_init_, press); const PVec sat = phaseSaturations(G, eqreg, cells, props, swat_init_, press);
@ -411,8 +416,8 @@ namespace Opm
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) { && props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid]; const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour]; const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]); const Vec rs = computeRs(G, cells, press[oilpos], temp, *(rs_func_[r]), sat[gaspos]);
const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]); const Vec rv = computeRs(G, cells, press[gaspos], temp, *(rv_func_[r]), sat[oilpos]);
copyFromRegion(rs, cells, rs_); copyFromRegion(rs, cells, rs_);
copyFromRegion(rv, cells, rv_); copyFromRegion(rv, cells, rv_);
} }

58
opm/core/simulator/initStateEquil_impl.hpp
View File

@ -106,11 +106,13 @@ namespace Opm
template <class Density> template <class Density>
class Water { class Water {
public: public:
Water(const Density& rho, Water(const double temp,
const Density& rho,
const int np, const int np,
const int ix, const int ix,
const double norm_grav) const double norm_grav)
: rho_(rho) : temp_(temp)
, rho_(rho)
, svol_(np, 0) , svol_(np, 0)
, ix_(ix) , ix_(ix)
, g_(norm_grav) , g_(norm_grav)
@ -126,6 +128,7 @@ namespace Opm
} }
private: private:
const double temp_;
const Density& rho_; const Density& rho_;
std::vector<double> svol_; std::vector<double> svol_;
const int ix_; const int ix_;
@ -134,7 +137,7 @@ namespace Opm
double double
density(const double press) const density(const double press) const
{ {
const std::vector<double>& rho = rho_(press, svol_); const std::vector<double>& rho = rho_(press, temp_, svol_);
return rho[ix_]; return rho[ix_];
} }
@ -143,13 +146,15 @@ namespace Opm
template <class Density, class RS> template <class Density, class RS>
class Oil { class Oil {
public: public:
Oil(const Density& rho, Oil(const double temp,
const Density& rho,
const RS& rs, const RS& rs,
const int np, const int np,
const int oix, const int oix,
const int gix, const int gix,
const double norm_grav) const double norm_grav)
: rho_(rho) : temp_(temp)
, rho_(rho)
, rs_(rs) , rs_(rs)
, svol_(np, 0) , svol_(np, 0)
, oix_(oix) , oix_(oix)
@ -167,6 +172,7 @@ namespace Opm
} }
private: private:
const double temp_;
const Density& rho_; const Density& rho_;
const RS& rs_; const RS& rs_;
mutable std::vector<double> svol_; mutable std::vector<double> svol_;
@ -179,10 +185,10 @@ namespace Opm
const double press) const const double press) const
{ {
if (gix_ >= 0) { if (gix_ >= 0) {
svol_[gix_] = rs_(depth, press); svol_[gix_] = rs_(depth, press, temp_);
} }
const std::vector<double>& rho = rho_(press, svol_); const std::vector<double>& rho = rho_(press, temp_, svol_);
return rho[oix_]; return rho[oix_];
} }
}; };
@ -190,13 +196,15 @@ namespace Opm
template <class Density, class RV> template <class Density, class RV>
class Gas { class Gas {
public: public:
Gas(const Density& rho, Gas(const double temp,
const Density& rho,
const RV& rv, const RV& rv,
const int np, const int np,
const int gix, const int gix,
const int oix, const int oix,
const double norm_grav) const double norm_grav)
: rho_(rho) : temp_(temp)
, rho_(rho)
, rv_(rv) , rv_(rv)
, svol_(np, 0) , svol_(np, 0)
, gix_(gix) , gix_(gix)
@ -214,6 +222,7 @@ namespace Opm
} }
private: private:
const double temp_;
const Density& rho_; const Density& rho_;
const RV& rv_; const RV& rv_;
mutable std::vector<double> svol_; mutable std::vector<double> svol_;
@ -226,10 +235,10 @@ namespace Opm
const double press) const const double press) const
{ {
if (oix_ >= 0) { if (oix_ >= 0) {
svol_[oix_] = rv_(depth, press); svol_[oix_] = rv_(depth, press, temp_);
} }
const std::vector<double>& rho = rho_(press, svol_); const std::vector<double>& rho = rho_(press, temp_, svol_);
return rho[gix_]; return rho[gix_];
} }
}; };
@ -333,7 +342,8 @@ namespace Opm
const PhaseUsage& pu = reg.phaseUsage(); const PhaseUsage& pu = reg.phaseUsage();
const int wix = PhaseIndex::water(pu); const int wix = PhaseIndex::water(pu);
ODE drho(reg.densityCalculator(), pu.num_phases, wix, grav); const double T = 273.15 + 20; // standard temperature for now
ODE drho(T, reg.densityCalculator(), pu.num_phases, wix, grav);
const double z0 = reg.zwoc(); const double z0 = reg.zwoc();
const double p0 = po_woc - reg.pcow_woc(); // Pcow = Po - Pw const double p0 = po_woc - reg.pcow_woc(); // Pcow = Po - Pw
@ -373,7 +383,9 @@ namespace Opm
const int oix = PhaseIndex::oil(pu); const int oix = PhaseIndex::oil(pu);
const int gix = PhaseIndex::gas(pu); const int gix = PhaseIndex::gas(pu);
ODE drho(reg.densityCalculator(), const double T = 273.15 + 20; // standard temperature for now
ODE drho(T,
reg.densityCalculator(),
reg.dissolutionCalculator(), reg.dissolutionCalculator(),
pu.num_phases, oix, gix, grav); pu.num_phases, oix, gix, grav);
@ -426,7 +438,9 @@ namespace Opm
const int gix = PhaseIndex::gas(pu); const int gix = PhaseIndex::gas(pu);
const int oix = PhaseIndex::oil(pu); const int oix = PhaseIndex::oil(pu);
ODE drho(reg.densityCalculator(), const double T = 273.15 + 20; // standard temperature for now
ODE drho(T,
reg.densityCalculator(),
reg.evaporationCalculator(), reg.evaporationCalculator(),
pu.num_phases, gix, oix, grav); pu.num_phases, gix, oix, grav);
@ -573,8 +587,16 @@ namespace Opm
return press; return press;
} }
template <class Region,
class CellRange>
std::vector<double>
temperature(const UnstructuredGrid& G,
const Region& reg,
const CellRange& cells)
{
// use the standard temperature for everything for now
return std::vector<double>(cells.size(), 273.15 + 20.0);
}
template <class Region, class CellRange> template <class Region, class CellRange>
std::vector< std::vector<double> > std::vector< std::vector<double> >
@ -716,6 +738,7 @@ namespace Opm
* \param[in] cells Range that spans the cells of the current * \param[in] cells Range that spans the cells of the current
* equilibration region. * equilibration region.
* \param[in] oil_pressure Oil pressure for each cell in range. * \param[in] oil_pressure Oil pressure for each cell in range.
* \param[in] temperature Temperature for each cell in range.
* \param[in] rs_func Rs as function of pressure and depth. * \param[in] rs_func Rs as function of pressure and depth.
* \return Rs values, one for each cell in the 'cells' range. * \return Rs values, one for each cell in the 'cells' range.
*/ */
@ -723,6 +746,7 @@ namespace Opm
std::vector<double> computeRs(const UnstructuredGrid& grid, std::vector<double> computeRs(const UnstructuredGrid& grid,
const CellRangeType& cells, const CellRangeType& cells,
const std::vector<double> oil_pressure, const std::vector<double> oil_pressure,
const std::vector<double>& temperature,
const Miscibility::RsFunction& rs_func, const Miscibility::RsFunction& rs_func,
const std::vector<double> gas_saturation) const std::vector<double> gas_saturation)
{ {
@ -731,7 +755,7 @@ namespace Opm
int count = 0; int count = 0;
for (auto it = cells.begin(); it != cells.end(); ++it, ++count) { for (auto it = cells.begin(); it != cells.end(); ++it, ++count) {
const double depth = grid.cell_centroids[3*(*it) + 2]; const double depth = grid.cell_centroids[3*(*it) + 2];
rs[count] = rs_func(depth, oil_pressure[count], gas_saturation[count]); rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]);
} }
return rs; return rs;
} }

19
opm/core/simulator/initState_impl.hpp
View File

@ -174,7 +174,7 @@ namespace Opm
{ {
const BlackoilPropertiesInterface& props_; const BlackoilPropertiesInterface& props_;
Density(const BlackoilPropertiesInterface& props) : props_(props) {} Density(const BlackoilPropertiesInterface& props) : props_(props) {}
double operator()(const double pressure, const int phase) double operator()(const double pressure, const double temperature, const int phase)
{ {
assert(props_.numPhases() == 2); assert(props_.numPhases() == 2);
const double surfvol[2][2] = { { 1.0, 0.0 }, const double surfvol[2][2] = { { 1.0, 0.0 },
@ -182,7 +182,7 @@ namespace Opm
// We do not handle multi-region PVT/EQUIL at this point. // We do not handle multi-region PVT/EQUIL at this point.
const int* cells = 0; const int* cells = 0;
double A[4] = { 0.0 }; double A[4] = { 0.0 };
props_.matrix(1, &pressure, surfvol[phase], cells, A, 0); props_.matrix(1, &pressure, &temperature, surfvol[phase], cells, A, 0);
double rho[2] = { 0.0 }; double rho[2] = { 0.0 };
props_.density(1, A, cells, rho); props_.density(1, A, cells, rho);
return rho[phase]; return rho[phase];
@ -233,11 +233,12 @@ namespace Opm
int phase = (datum_z > woc) ? 0 : 1; int phase = (datum_z > woc) ? 0 : 1;
int num_steps = int(std::ceil(std::fabs(woc - datum_z)/hmin)); int num_steps = int(std::ceil(std::fabs(woc - datum_z)/hmin));
double pval = datum_p; double pval = datum_p;
double Tval = 273.15 + 20; // standard temperature
double zval = datum_z; double zval = datum_z;
double h = (woc - datum_z)/double(num_steps); double h = (woc - datum_z)/double(num_steps);
for (int i = 0; i < num_steps; ++i) { for (int i = 0; i < num_steps; ++i) {
zval += h; zval += h;
const double dp = rho(pval, phase)*gravity; const double dp = rho(pval, Tval, phase)*gravity;
pval += h*dp; pval += h*dp;
press_by_z[zval] = pval; press_by_z[zval] = pval;
} }
@ -251,7 +252,7 @@ namespace Opm
h = (z_end - datum_z)/double(num_steps); h = (z_end - datum_z)/double(num_steps);
for (int i = 0; i < num_steps; ++i) { for (int i = 0; i < num_steps; ++i) {
zval += h; zval += h;
const double dp = rho(pval, phase)*gravity; const double dp = rho(pval, Tval, phase)*gravity;
pval += h*dp; pval += h*dp;
press_by_z[zval] = pval; press_by_z[zval] = pval;
} }
@ -265,7 +266,7 @@ namespace Opm
h = (z_end - datum_z)/double(num_steps); h = (z_end - datum_z)/double(num_steps);
for (int i = 0; i < num_steps; ++i) { for (int i = 0; i < num_steps; ++i) {
zval += h; zval += h;
const double dp = rho(pval, phase)*gravity; const double dp = rho(pval, Tval, phase)*gravity;
pval += h*dp; pval += h*dp;
press_by_z[zval] = pval; press_by_z[zval] = pval;
} }
@ -699,7 +700,7 @@ namespace Opm
// Assuming that using the saturation as z argument here does not change // Assuming that using the saturation as z argument here does not change
// the outcome. This is not guaranteed unless we have only a single phase // the outcome. This is not guaranteed unless we have only a single phase
// per cell. // per cell.
props.matrix(nc, &state.pressure()[0], &state.surfacevol()[0], &allcells[0], &allA[0], 0); props.matrix(nc, &state.pressure()[0], &state.temperature()[0], &state.surfacevol()[0], &allcells[0], &allA[0], 0);
for (int c = 0; c < nc; ++c) { for (int c = 0; c < nc; ++c) {
// Using z = As // Using z = As
double* z = &state.surfacevol()[c*np]; double* z = &state.surfacevol()[c*np];
@ -783,7 +784,7 @@ namespace Opm
z_init[c*np + p] = z_tmp; z_init[c*np + p] = z_tmp;
} }
} }
props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_a[0], 0); props.matrix(number_of_cells, &state.pressure()[0], &state.temperature()[0], &z_init[0], &allcells[0], &allA_a[0], 0);
// Liquid phase // Liquid phase
if(pu.phase_used[BlackoilPhases::Liquid]){ if(pu.phase_used[BlackoilPhases::Liquid]){
@ -805,7 +806,7 @@ namespace Opm
} }
} }
} }
props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_l[0], 0); props.matrix(number_of_cells, &state.pressure()[0], &state.temperature()[0], &z_init[0], &allcells[0], &allA_l[0], 0);
if(pu.phase_used[BlackoilPhases::Vapour]){ if(pu.phase_used[BlackoilPhases::Vapour]){
for (int c = 0; c < number_of_cells ; ++c){ for (int c = 0; c < number_of_cells ; ++c){
@ -826,7 +827,7 @@ namespace Opm
} }
} }
} }
props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_v[0], 0); props.matrix(number_of_cells, &state.pressure()[0], &state.temperature()[0], &z_init[0], &allcells[0], &allA_v[0], 0);
for (int c = 0; c < number_of_cells; ++c) { for (int c = 0; c < number_of_cells; ++c) {
// Using z = As // Using z = As

5
opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp
View File

@ -80,6 +80,7 @@ namespace Opm
void TransportSolverCompressibleTwophaseReorder::solve(const double* darcyflux, void TransportSolverCompressibleTwophaseReorder::solve(const double* darcyflux,
const double* pressure, const double* pressure,
const double* temperature,
const double* porevolume0, const double* porevolume0,
const double* porevolume, const double* porevolume,
const double* source, const double* source,
@ -95,8 +96,8 @@ namespace Opm
dt_ = dt; dt_ = dt;
toWaterSat(saturation, saturation_); toWaterSat(saturation, saturation_);
props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); props_.viscosity(props_.numCells(), pressure, temperature, NULL, &allcells_[0], &visc_[0], NULL);
props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL); props_.matrix(props_.numCells(), pressure, temperature, NULL, &allcells_[0], &A_[0], NULL);
// Check immiscibility requirement (only done for first cell). // Check immiscibility requirement (only done for first cell).
if (A_[1] != 0.0 || A_[2] != 0.0) { if (A_[1] != 0.0 || A_[2] != 0.0) {

2
opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp
View File

@ -48,6 +48,7 @@ namespace Opm
/// Solve for saturation at next timestep. /// Solve for saturation at next timestep.
/// \param[in] darcyflux Array of signed face fluxes. /// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] pressure Array of cell pressures /// \param[in] pressure Array of cell pressures
/// \param[in] temperature Array of cell temperatures
/// \param[in] surfacevol0 Array of surface volumes at start of timestep /// \param[in] surfacevol0 Array of surface volumes at start of timestep
/// \param[in] porevolume0 Array of pore volumes at start of timestep. /// \param[in] porevolume0 Array of pore volumes at start of timestep.
/// \param[in] porevolume Array of pore volumes at end of timestep. /// \param[in] porevolume Array of pore volumes at end of timestep.
@ -57,6 +58,7 @@ namespace Opm
/// \param[in, out] surfacevol Surface volume densities for each phase. /// \param[in, out] surfacevol Surface volume densities for each phase.
void solve(const double* darcyflux, void solve(const double* darcyflux,
const double* pressure, const double* pressure,
const double* temperature,
const double* porevolume0, const double* porevolume0,
const double* porevolume, const double* porevolume,
const double* source, const double* source,

2
opm/core/utility/miscUtilities.cpp
View File

@ -636,7 +636,7 @@ namespace Opm
double mob[max_np]; double mob[max_np];
props.relperm(1, &s[np*cell], &cell, mob, 0); props.relperm(1, &s[np*cell], &cell, mob, 0);
double visc[max_np]; double visc[max_np];
props.viscosity(1, &p[cell], &z[np*cell], &cell, visc, 0); props.viscosity(1, &p[cell], 0, &z[np*cell], &cell, visc, 0);
double tmob = 0; double tmob = 0;
for(int i = 0; i < np; ++i) { for(int i = 0; i < np; ++i) {
mob[i] /= visc[i]; mob[i] /= visc[i];

21
opm/core/utility/miscUtilitiesBlackoil.cpp
View File

@ -74,6 +74,7 @@ namespace Opm
OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells."); OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells.");
} }
const std::vector<double>& press = state.pressure(); const std::vector<double>& press = state.pressure();
const std::vector<double>& temp = state.temperature();
const std::vector<double>& s = state.saturation(); const std::vector<double>& s = state.saturation();
const std::vector<double>& z = state.surfacevol(); const std::vector<double>& z = state.surfacevol();
std::fill(injected, injected + np, 0.0); std::fill(injected, injected + np, 0.0);
@ -94,8 +95,8 @@ namespace Opm
const double flux = -transport_src[c]*dt; const double flux = -transport_src[c]*dt;
const double* sat = &s[np*c]; const double* sat = &s[np*c];
props.relperm(1, sat, &c, &mob[0], 0); props.relperm(1, sat, &c, &mob[0], 0);
props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0); props.viscosity(1, &press[c], &temp[c], &z[np*c], &c, &visc[0], 0);
props.matrix(1, &press[c], &z[np*c], &c, &A[0], 0); props.matrix(1, &press[c], &temp[c], &z[np*c], &c, &A[0], 0);
double totmob = 0.0; double totmob = 0.0;
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
mob[p] /= visc[p]; mob[p] /= visc[p];
@ -121,19 +122,21 @@ namespace Opm
/// @param[in] props rock and fluid properties /// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated /// @param[in] cells cells with which the saturation values are associated
/// @param[in] p pressure (one value per cell) /// @param[in] p pressure (one value per cell)
/// @param[in] temp temperature (one value per cell)
/// @param[in] z surface-volume values (for all P phases) /// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all phases) /// @param[in] s saturation values (for all phases)
/// @param[out] totmob total mobilities. /// @param[out] totmob total mobilities.
void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props, void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells, const std::vector<int>& cells,
const std::vector<double>& press, const std::vector<double>& press,
const std::vector<double>& temp,
const std::vector<double>& z, const std::vector<double>& z,
const std::vector<double>& s, const std::vector<double>& s,
std::vector<double>& totmob) std::vector<double>& totmob)
{ {
std::vector<double> pmobc; std::vector<double> pmobc;
computePhaseMobilities(props, cells, press, z, s, pmobc); computePhaseMobilities(props, cells, press, temp, z, s, pmobc);
const std::size_t np = props.numPhases(); const std::size_t np = props.numPhases();
const std::vector<int>::size_type nc = cells.size(); const std::vector<int>::size_type nc = cells.size();
@ -193,12 +196,14 @@ namespace Opm
/// @param[in] props rock and fluid properties /// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated /// @param[in] cells cells with which the saturation values are associated
/// @param[in] p pressure (one value per cell) /// @param[in] p pressure (one value per cell)
/// @param[in] T temperature (one value per cell)
/// @param[in] z surface-volume values (for all P phases) /// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all phases) /// @param[in] s saturation values (for all phases)
/// @param[out] pmobc phase mobilities (for all phases). /// @param[out] pmobc phase mobilities (for all phases).
void computePhaseMobilities(const Opm::BlackoilPropertiesInterface& props, void computePhaseMobilities(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells, const std::vector<int>& cells,
const std::vector<double>& p, const std::vector<double>& p,
const std::vector<double>& T,
const std::vector<double>& z, const std::vector<double>& z,
const std::vector<double>& s, const std::vector<double>& s,
std::vector<double>& pmobc) std::vector<double>& pmobc)
@ -209,7 +214,7 @@ namespace Opm
assert(int(s.size()) == nc * np); assert(int(s.size()) == nc * np);
std::vector<double> mu(nc*np); std::vector<double> mu(nc*np);
props.viscosity(nc, &p[0], &z[0], &cells[0], &mu[0], 0); props.viscosity(nc, &p[0], &T[0], &z[0], &cells[0], &mu[0], 0);
pmobc.clear(); pmobc.clear();
pmobc.resize(nc*np, 0.0); pmobc.resize(nc*np, 0.0);
@ -227,19 +232,21 @@ namespace Opm
/// @param[in] props rock and fluid properties /// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated /// @param[in] cells cells with which the saturation values are associated
/// @param[in] p pressure (one value per cell) /// @param[in] p pressure (one value per cell)
/// @param[in] T temperature (one value per cell)
/// @param[in] z surface-volume values (for all P phases) /// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all phases) /// @param[in] s saturation values (for all phases)
/// @param[out] fractional_flow the fractional flow for each phase for each cell. /// @param[out] fractional_flow the fractional flow for each phase for each cell.
void computeFractionalFlow(const Opm::BlackoilPropertiesInterface& props, void computeFractionalFlow(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells, const std::vector<int>& cells,
const std::vector<double>& p, const std::vector<double>& p,
const std::vector<double>& T,
const std::vector<double>& z, const std::vector<double>& z,
const std::vector<double>& s, const std::vector<double>& s,
std::vector<double>& fractional_flows) std::vector<double>& fractional_flows)
{ {
const int num_phases = props.numPhases(); const int num_phases = props.numPhases();
computePhaseMobilities(props, cells, p, z, s, fractional_flows); computePhaseMobilities(props, cells, p, T, z, s, fractional_flows);
for (std::vector<int>::size_type i = 0; i < cells.size(); ++i) { for (std::vector<int>::size_type i = 0; i < cells.size(); ++i) {
double phase_sum = 0.0; double phase_sum = 0.0;
@ -299,7 +306,7 @@ namespace Opm
//std::vector<double> res_vol(np); //std::vector<double> res_vol(np);
const std::vector<double>& z = state.surfacevol(); const std::vector<double>& z = state.surfacevol();
props.matrix(nc, &state.pressure()[0], &z[0], &allcells[0], &allA[0], 0); props.matrix(nc, &state.pressure()[0], &state.temperature()[0], &z[0], &allcells[0], &allA[0], 0);
// Linear solver. // Linear solver.
MAT_SIZE_T n = np; MAT_SIZE_T n = np;
@ -388,7 +395,7 @@ namespace Opm
} else { } else {
assert(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6); assert(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
perf_rate *= comp_frac[0]; // Water reservoir volume rate. perf_rate *= comp_frac[0]; // Water reservoir volume rate.
props.matrix(1, &well_state.perfPress()[perf], comp_frac, &perf_cell, &A[0], 0); props.matrix(1, &well_state.perfPress()[perf], &well_state.temperature()[w], comp_frac, &perf_cell, &A[0], 0);
perf_rate *= A[0]; // Water surface volume rate. perf_rate *= A[0]; // Water surface volume rate.
} }
} }

4
opm/core/utility/miscUtilitiesBlackoil.hpp
View File

@ -92,12 +92,14 @@ namespace Opm
/// @param[in] props rock and fluid properties /// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated /// @param[in] cells cells with which the saturation values are associated
/// @param[in] p pressure (one value per cell) /// @param[in] p pressure (one value per cell)
/// @param[in] T temperature (one value per cell)
/// @param[in] z surface-volume values (for all P phases) /// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all phases) /// @param[in] s saturation values (for all phases)
/// @param[out] pmobc phase mobilities (for all phases). /// @param[out] pmobc phase mobilities (for all phases).
void computePhaseMobilities(const Opm::BlackoilPropertiesInterface& props, void computePhaseMobilities(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells, const std::vector<int>& cells,
const std::vector<double>& p, const std::vector<double>& p,
const std::vector<double>& T,
const std::vector<double>& z, const std::vector<double>& z,
const std::vector<double>& s, const std::vector<double>& s,
std::vector<double>& pmobc); std::vector<double>& pmobc);
@ -107,12 +109,14 @@ namespace Opm
/// @param[in] props rock and fluid properties /// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated /// @param[in] cells cells with which the saturation values are associated
/// @param[in] p pressure (one value per cell) /// @param[in] p pressure (one value per cell)
/// @param[in] T temperature(one value per cell)
/// @param[in] z surface-volume values (for all P phases) /// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all phases) /// @param[in] s saturation values (for all phases)
/// @param[out] fractional_flow the fractional flow for each phase for each cell. /// @param[out] fractional_flow the fractional flow for each phase for each cell.
void computeFractionalFlow(const Opm::BlackoilPropertiesInterface& props, void computeFractionalFlow(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells, const std::vector<int>& cells,
const std::vector<double>& p, const std::vector<double>& p,
const std::vector<double>& T,
const std::vector<double>& z, const std::vector<double>& z,
const std::vector<double>& s, const std::vector<double>& s,
std::vector<double>& fractional_flows); std::vector<double>& fractional_flows);

22
tests/test_blackoilfluid.cpp
View File

@ -102,7 +102,7 @@ std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, Opm
return props_; return props_;
} }
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, void testmu(const double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<double> T, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<PvtInterface> > 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> mu(n);
@ -116,8 +116,8 @@ void testmu(const double reltol, int n, int np, const std::vector<int> &pvtTable
// test mu // test mu
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
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], &T[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &z[0], &mu[0]); props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &T[0], &z[0], &mu[0]);
dmudp_diff = (mu_new[1]-mu_new[0])/(p[1]-p[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]); 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]); dmudp_diff_u = (mu_new[4]-mu_new[3])/(p[4]-p[3]);
@ -138,7 +138,7 @@ void testmu(const double reltol, int n, int np, const std::vector<int> &pvtTable
} }
} }
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, void testb(const double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<double> T, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition) std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{ {
// test b // test b
@ -155,8 +155,8 @@ void testb(const double reltol, int n, int np, const std::vector<int> &pvtTableI
double dbdr_diff_u; double dbdr_diff_u;
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
props_[phase]->b(n, &pvtTableIdx[0], &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]); props_[phase]->b(n, &pvtTableIdx[0], &p[0], &T[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]); props_[phase]->dBdp(n, &pvtTableIdx[0], &p[0], &T[0], &z[0], &B[0], &dBdp[0]);
dbdp_diff = (b[1]-b[0])/(p[1]-p[0]); dbdp_diff = (b[1]-b[0])/(p[1]-p[0]);
dbdr_diff = (b[2]-b[0])/(r[2]-r[0]); dbdr_diff = (b[2]-b[0])/(r[2]-r[0]);
dbdp_diff_u = (b[4]-b[3])/(p[4]-p[3]); dbdp_diff_u = (b[4]-b[3])/(p[4]-p[3]);
@ -253,6 +253,7 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
const double reltolpermu = 1e-1; const double reltolpermu = 1e-1;
std::vector<double> p(n); std::vector<double> p(n);
std::vector<double> T(n, 273.15 + 20);
std::vector<double> r(n); std::vector<double> r(n);
std::vector<PhasePresence> condition(n); std::vector<PhasePresence> condition(n);
std::vector<double> z(n * np); std::vector<double> z(n * np);
@ -292,9 +293,9 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
} }
testmu(reltolpermu, n, np, pvtRegionIdx, p, r,z, props_, condition); testmu(reltolpermu, n, np, pvtRegionIdx, p, T, r,z, props_, condition);
testb(reltolper,n,np,pvtRegionIdx,p,r,z,props_,condition); testb(reltolper,n,np,pvtRegionIdx,p,T,r,z,props_,condition);
testrsSat(reltolper,n,np,pvtRegionIdx,p,props_); testrsSat(reltolper,n,np,pvtRegionIdx,p,props_);
@ -332,6 +333,7 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
const double reltolpermu = 1e-1; const double reltolpermu = 1e-1;
std::vector<double> p(n); std::vector<double> p(n);
std::vector<double> T(n, 273.15+20);
std::vector<double> r(n); std::vector<double> r(n);
std::vector<PhasePresence> condition(n); std::vector<PhasePresence> condition(n);
std::vector<double> z(n * np); std::vector<double> z(n * np);
@ -371,9 +373,9 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
} }
testmu(reltolpermu, n, np, pvtRegionIdx, p, r,z, props_, condition); testmu(reltolpermu, n, np, pvtRegionIdx, p,T, r,z, props_, condition);
testb(reltolper,n,np,pvtRegionIdx,p,r,z,props_,condition); testb(reltolper,n,np,pvtRegionIdx,p,T,r,z,props_,condition);
testrsSat(reltolper,n,np,pvtRegionIdx,p,props_); testrsSat(reltolper,n,np,pvtRegionIdx,p,props_);