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

PVT properties: allow them to be temperature dependent

Browse Source 
Note that this patch does not introduce any real temperature
dependence but only changes the APIs for the viscosity and for the
density related methods. Note that I also don't like the fact that
this requires so many changes to so many files, but with the current
design of the property classes I cannot see a way to avoid this...
This commit is contained in:
Andreas Lauser 2014-11-20 12:15:01 +01:00
parent 960dac3caf
commit e530ce03d8
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) {
const int cell = wells_->well_cells[j];
const double cell_depth = grid_.cell_centroids[dim * cell + dim - 1];
props_.matrix(1, &state.pressure()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0);
props_.matrix(1, &state.pressure()[cell], &state.temperature()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0);
props_.density(1, &A[0], &cell, &rho[0]);
for (int phase = 0; phase < np; ++phase) {
const double s_phase = state.saturation()[np*cell + phase];
@ -309,13 +309,14 @@ namespace Opm
const int nc = grid_.number_of_cells;
const int np = props_.numPhases();
const double* cell_p = &state.pressure()[0];
const double* cell_T = &state.temperature()[0];
const double* cell_z = &state.surfacevol()[0];
const double* cell_s = &state.saturation()[0];
cell_A_.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);
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);
props_.relperm(nc, cell_s, &allcells_[0], &cell_phasemob_[0], 0);
std::transform(cell_phasemob_.begin(), cell_phasemob_.end(),
@ -481,13 +482,14 @@ namespace Opm
} else {
const double bhp = well_state.bhp()[w];
double perf_p = bhp + wellperf_wdp_[j];
const double perf_T = well_state.temperature()[w];
// Hack warning: comp_frac is used as a component
// surface-volume variable in calls to matrix() and
// viscosity(), but as a saturation in the call to
// relperm(). This is probably ok as long as injectors
// only inject pure fluids.
props_.matrix(1, &perf_p, comp_frac, &c, wpA, NULL);
props_.viscosity(1, &perf_p, comp_frac, &c, &mu[0], NULL);
props_.matrix(1, &perf_p, &perf_T, comp_frac, &c, wpA, 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);
props_.relperm (1, comp_frac, &c, wpM , NULL);
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] 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] 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.
@ -98,6 +99,7 @@ namespace Opm
/// array must be valid before calling.
void BlackoilPropertiesBasic::viscosity(const int n,
const double* p,
const double* T,
const double* z,
const int* /*cells*/,
double* mu,
@ -106,12 +108,13 @@ namespace Opm
if (dmudp) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesBasic::viscosity() -- derivatives of viscosity not yet implemented.");
} else {
pvt_.mu(n, p, z, mu);
pvt_.mu(n, p, T, z, mu);
}
}
/// \param[in] n Number of data points.
/// \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] 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.
@ -121,7 +124,8 @@ namespace Opm
/// array must be valid before calling. The matrices are output
/// in Fortran order.
void BlackoilPropertiesBasic::matrix(const int n,
const double* /*p*/,
const double* p,
const double* T,
const double* /*z*/,
const int* /*cells*/,
double* A,
@ -130,7 +134,7 @@ namespace Opm
const int np = numPhases();
assert(np <= 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
// #pragma omp parallel for
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] 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] 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.
@ -90,6 +91,7 @@ namespace Opm
/// array must be valid before calling.
virtual void viscosity(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* mu,
@ -97,6 +99,7 @@ namespace Opm
/// \param[in] n Number of data points.
/// \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] 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.
@ -107,6 +110,7 @@ namespace Opm
/// in Fortran order.
virtual void matrix(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
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] 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] 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.
@ -97,6 +98,7 @@ namespace Opm
/// array must be valid before calling.
void BlackoilPropertiesFromDeck::viscosity(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* mu,
@ -111,12 +113,13 @@ namespace Opm
for (int i = 0; i < n; ++ 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] 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] 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.
@ -127,6 +130,7 @@ namespace Opm
/// in Fortran order.
void BlackoilPropertiesFromDeck::matrix(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* A,
@ -144,10 +148,10 @@ namespace Opm
if (dAdp) {
dB_.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]);
} 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]);
}
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] 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] 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.
@ -132,6 +133,7 @@ namespace Opm
/// array must be valid before calling.
virtual void viscosity(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* mu,
@ -139,6 +141,7 @@ namespace Opm
/// \param[in] n Number of data points.
/// \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] 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.
@ -149,6 +152,7 @@ namespace Opm
/// in Fortran order.
virtual void matrix(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
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] 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] 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.
@ -77,6 +78,7 @@ namespace Opm
/// array must be valid before calling.
virtual void viscosity(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* mu,
@ -84,6 +86,7 @@ namespace Opm
/// \param[in] n Number of data points.
/// \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] 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.
@ -94,6 +97,7 @@ namespace Opm
/// in Fortran order.
virtual void matrix(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
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() << ").");
}
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,
@ -64,7 +64,7 @@ namespace Opm
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]);
pvt_.mu(1, 0, 0, 0, &viscosity_[0]);
}
IncompPropertiesBasic::~IncompPropertiesBasic()

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

@ -161,12 +161,13 @@ namespace Opm
void BlackoilPvtProperties::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const
{
data1_.resize(n);
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
for (int i = 0; i < n; ++i) {
output_mu[phase_usage_.num_phases*i + phase] = data1_[i];
@ -177,12 +178,13 @@ namespace Opm
void BlackoilPvtProperties::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const
{
data1_.resize(n);
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
for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i];
@ -193,6 +195,7 @@ namespace Opm
void BlackoilPvtProperties::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const
@ -200,7 +203,7 @@ namespace Opm
data1_.resize(n);
data2_.resize(n);
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
for (int i = 0; i < n; ++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().
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,
const int *pvtTableIdx,
const double* p,
const double* T,
const double* z,
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,
const int *pvtTableIdx,
const double* p,
const double* T,
const double* z,
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,
const int *pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
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,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*z*/,
double* output_mu) const
{
@ -124,6 +125,7 @@ namespace Opm
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
double* output_mu,
double* output_dmudp,
@ -144,6 +146,7 @@ namespace Opm
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_mu,
@ -165,6 +168,7 @@ namespace Opm
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*z*/,
double* output_B) const
{
@ -180,6 +184,7 @@ namespace Opm
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
@ -197,6 +202,7 @@ namespace Opm
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
double* output_b,
double* output_dbdp,
@ -220,6 +226,7 @@ namespace Opm
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_b,

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

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

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

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

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

@ -49,64 +49,71 @@ namespace Opm
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,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
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).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
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'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
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,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
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,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
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).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
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'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
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.
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)
/// or pressure (p) and gas resolution factor (r).
/// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z)
/// or pressure (p), temperature (T) and gas resolution factor (r).
/// For all the virtual methods, the following apply:
/// - pvtRegionIdx is an array of size n and represents the
/// index of the PVT table which should be used to calculate
@ -54,38 +54,42 @@ namespace Opm
/// - Output arrays shall be of size n, and must be valid before
/// 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,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
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).
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
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'.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
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,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const = 0;
@ -93,25 +97,28 @@ namespace Opm
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
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).
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
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'.
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_b,

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

@ -104,6 +104,7 @@ namespace Opm
void PvtLiveGas::mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* z,
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*/,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*p*/,
const double* /*T*/,
const double* /*r*/,
double* /*output_mu*/,
double* /*output_dmudp*/,
@ -128,10 +130,11 @@ namespace Opm
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,
const int* pvtRegionIdx,
const double* p,
const double* p,
const double* /*T*/,
const double* r,
const PhasePresence* cond,
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,
const int* pvtRegionIdx,
const double* p,
const double* p,
const double* /*T*/,
const double* z,
double* output_B) const
{
@ -181,8 +185,9 @@ namespace Opm
/// Formation volume factor and p-derivative as functions of p and z.
void PvtLiveGas::dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* z,
double* output_B,
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*/,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*p*/,
const double* /*T*/,
const double* /*r*/,
double* /*output_b*/,
double* /*output_dbdp*/,
@ -206,10 +212,11 @@ namespace Opm
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,
const int* pvtRegionIdx,
const double* p,
const double* p,
const double* /*T*/,
const double* r,
const PhasePresence* cond,
double* output_b,

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

@ -29,8 +29,8 @@
namespace Opm
{
/// 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)
/// or pressure (p) and gas resolution factor (r).
/// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z)
/// or pressure (p), temperature (T) and gas resolution factor (r).
/// For all the virtual methods, the following apply: p, r and z
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
@ -41,64 +41,71 @@ namespace Opm
PvtLiveGas(const std::vector<Opm::PvtgTable>& pvtgTables);
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,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
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).
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
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'.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
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,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
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,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
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).
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
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'.
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
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,
const int* pvtRegionIdx,
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,
const int* pvtTableIdx,
const double* p,
const double* p,
const double* /*T*/,
const double* z,
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,
const int* pvtTableIdx,
const double* p,
const double* p,
const double* /*T*/,
const double* r,
double* output_mu,
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,
const int* pvtTableIdx,
const double* p,
const double* p,
const double* /*T*/,
const double* r,
const PhasePresence* cond,
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,
const int* pvtTableIdx,
const double* p,
const double* p,
const double* /*T*/,
const double* z,
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,
const int* pvtTableIdx,
const double* p,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* z,
double* output_B,
double* output_dBdp) const
@ -255,8 +260,9 @@ namespace Opm
}
void PvtLiveOil::b(const int n,
const int* pvtTableIdx,
const double* p,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* r,
double* output_b,
double* output_dbdp,
@ -275,8 +281,9 @@ namespace Opm
}
void PvtLiveOil::b(const int n,
const int* pvtTableIdx,
const double* p,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* r,
const PhasePresence* cond,
double* output_b,

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

@ -30,8 +30,8 @@
namespace Opm
{
/// 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)
/// or pressure (p) and gas resolution factor (r).
/// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z)
/// or pressure (p), temperature (T) and gas resolution factor (r).
/// For all the virtual methods, the following apply: p, r and z
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
@ -42,64 +42,71 @@ namespace Opm
PvtLiveOil(const std::vector<Opm::PvtoTable>& pvtoTables);
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,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
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).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
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'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
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,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
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,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
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).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
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'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_b,

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

@ -113,6 +113,7 @@ namespace Opm
void PvtPropertiesBasic::mu(const int n,
const double* /*p*/,
const double* /*T*/,
const double* /*z*/,
double* output_mu) const
{
@ -127,6 +128,7 @@ namespace Opm
void PvtPropertiesBasic::B(const int n,
const double* /*p*/,
const double* /*T*/,
const double* /*z*/,
double* output_B) const
{
@ -141,6 +143,7 @@ namespace Opm
void PvtPropertiesBasic::dBdp(const int n,
const double* /*p*/,
const double* /*T*/,
const double* /*z*/,
double* output_B,
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.
/// 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
/// before calling the method.
/// NOTE: This class is intentionally similar to BlackoilPvtProperties.
@ -62,21 +62,24 @@ namespace Opm
/// \return Array of size numPhases().
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,
const double* p,
const double* T,
const double* z,
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,
const double* p,
const double* T,
const double* z,
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,
const double* p,
const double* T,
const double* z,
double* output_B,
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] T Temperature.
*
* \param[in] z Surface volumes of all phases.
*
* \return Phase densities at phase point.
*/
std::vector<double>
operator()(const double p,
const double T,
const std::vector<double>& z) const
{
const int np = props_.numPhases();
@ -136,7 +139,7 @@ namespace Opm
assert (z.size() == std::vector<double>::size_type(np));
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);
props_.density(1, &A[0], &c_[0], &rho[0]);
@ -171,11 +174,15 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
virtual double operator()(const double depth,
const double press,
const double temp,
const double sat = 0.0) const = 0;
};
@ -194,6 +201,9 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. In "no mixing
* policy", this is identically zero.
@ -201,6 +211,7 @@ namespace Opm
double
operator()(const double /* depth */,
const double /* press */,
const double /* temp */,
const double /* sat */ = 0.0) const
{
return 0.0;
@ -247,18 +258,22 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double depth,
const double press,
const double temp,
const double sat_gas = 0.0) const
{
if (sat_gas > 0.0) {
return satRs(press);
return satRs(press, temp);
} 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];
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_.
// 1/Bo is in the oil row and column.
// 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
* value.
*
* \param[in] temp Temperature at which to calculate RV
* value.
*
* \return Vaporized oil-gas ratio (RV) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double depth,
const double press,
const double temp,
const double sat_oil = 0.0 ) const
{
if (sat_oil > 0.0) {
return satRv(press);
return satRv(press, temp);
} 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];
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_.
// 1/Bg is in the gas row and column.
// Recall also that it is stored in column-major order.
@ -382,15 +401,16 @@ namespace Opm
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \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)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
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
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double /* depth */,
const double press,
const double temp,
const double sat_gas = 0.0) const
{
if (sat_gas > 0.0) {
return satRs(press);
return satRs(press, temp);
} 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_;
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_.
// 1/Bo is in the oil row and column.
// Recall also that it is stored in column-major order.
@ -460,15 +484,16 @@ namespace Opm
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \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)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
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
* value.
*
* \param[in] temp Temperature at which to calculate RV
* value.
*
* \return Dissolved oil-gas ratio (RV) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double /*depth*/,
const double press,
const double temp,
const double sat_oil = 0.0) const
{
if (sat_oil > 0.0) {
return satRv(press);
return satRv(press, temp);
} 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_;
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_.
// 1/Bg is in the gas row and column.
// 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.
if (check_well_controls_) {
computeFractionalFlow(props_, allcells_,
state.pressure(), state.surfacevol(), state.saturation(),
state.pressure(), state.temperature(), state.surfacevol(), state.saturation(),
fractional_flows);
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
}
@ -445,7 +445,7 @@ namespace Opm
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
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,
state.saturation(), state.surfacevol());
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
equal = equal && vectorApproxEqual( pressure() , other.pressure() , epsilon);
equal = equal && vectorApproxEqual( temperature() , other.temperature() , epsilon);
equal = equal && vectorApproxEqual( facepressure() , other.facepressure() , epsilon);
equal = equal && vectorApproxEqual( faceflux() , other.faceflux() , 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;
press_.resize(number_of_cells, 0.0);
temp_.resize(number_of_cells, 273.15 + 20);
fpress_.resize(number_of_faces, 0.0);
flux_.resize(number_of_faces, 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_; }
std::vector<double>& pressure () { return press_ ; }
std::vector<double>& temperature () { return temp_ ; }
std::vector<double>& facepressure() { return fpress_; }
std::vector<double>& faceflux () { return flux_ ; }
std::vector<double>& saturation () { return sat_ ; }
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>& faceflux () const { return flux_ ; }
const std::vector<double>& saturation () const { return sat_ ; }
@ -56,6 +58,7 @@ namespace Opm
private:
int num_phases_;
std::vector<double> press_ ;
std::vector<double> temp_ ;
std::vector<double> fpress_;
std::vector<double> flux_ ;
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 np = wells->number_of_phases;
bhp_.resize(nw);
temperature_.resize(nw, 273.15 + 20); // standard temperature for now
wellrates_.resize(nw * np, 0.0);
for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
@ -110,6 +111,10 @@ namespace Opm
std::vector<double>& bhp() { 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.
std::vector<double>& wellRates() { return wellrates_; }
const std::vector<double>& wellRates() const { return wellrates_; }
@ -124,6 +129,7 @@ namespace Opm
private:
std::vector<double> bhp_;
std::vector<double> temperature_;
std::vector<double> wellrates_;
std::vector<double> perfrates_;
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
* equilibration region.
* \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.
* \return Rs values, one for each cell in the 'cells' range.
*/
@ -178,6 +179,7 @@ namespace Opm
std::vector<double> computeRs(const UnstructuredGrid& grid,
const CellRangeType& cells,
const std::vector<double> oil_pressure,
const std::vector<double>& temperature,
const Miscibility::RsFunction& rs_func,
const std::vector<double> gas_saturation);
@ -298,7 +300,8 @@ namespace Opm
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
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 {
@ -329,7 +332,8 @@ namespace Opm
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
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 {
@ -399,6 +403,7 @@ namespace Opm
props.phaseUsage());
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);
@ -411,8 +416,8 @@ namespace Opm
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
const Vec rs = computeRs(G, cells, press[oilpos], temp, *(rs_func_[r]), sat[gaspos]);
const Vec rv = computeRs(G, cells, press[gaspos], temp, *(rv_func_[r]), sat[oilpos]);
copyFromRegion(rs, cells, rs_);
copyFromRegion(rv, cells, rv_);
}

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

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

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

@ -174,7 +174,7 @@ namespace Opm
{
const BlackoilPropertiesInterface& 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);
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.
const int* cells = 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 };
props_.density(1, A, cells, rho);
return rho[phase];
@ -233,11 +233,12 @@ namespace Opm
int phase = (datum_z > woc) ? 0 : 1;
int num_steps = int(std::ceil(std::fabs(woc - datum_z)/hmin));
double pval = datum_p;
double Tval = 273.15 + 20; // standard temperature
double zval = datum_z;
double h = (woc - datum_z)/double(num_steps);
for (int i = 0; i < num_steps; ++i) {
zval += h;
const double dp = rho(pval, phase)*gravity;
const double dp = rho(pval, Tval, phase)*gravity;
pval += h*dp;
press_by_z[zval] = pval;
}
@ -251,7 +252,7 @@ namespace Opm
h = (z_end - datum_z)/double(num_steps);
for (int i = 0; i < num_steps; ++i) {
zval += h;
const double dp = rho(pval, phase)*gravity;
const double dp = rho(pval, Tval, phase)*gravity;
pval += h*dp;
press_by_z[zval] = pval;
}
@ -265,7 +266,7 @@ namespace Opm
h = (z_end - datum_z)/double(num_steps);
for (int i = 0; i < num_steps; ++i) {
zval += h;
const double dp = rho(pval, phase)*gravity;
const double dp = rho(pval, Tval, phase)*gravity;
pval += h*dp;
press_by_z[zval] = pval;
}
@ -699,7 +700,7 @@ namespace Opm
// 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
// 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) {
// Using z = As
double* z = &state.surfacevol()[c*np];
@ -783,7 +784,7 @@ namespace Opm
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
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]){
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) {
// Using z = As

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

@ -80,6 +80,7 @@ namespace Opm
void TransportSolverCompressibleTwophaseReorder::solve(const double* darcyflux,
const double* pressure,
const double* temperature,
const double* porevolume0,
const double* porevolume,
const double* source,
@ -95,8 +96,8 @@ namespace Opm
dt_ = dt;
toWaterSat(saturation, saturation_);
props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL);
props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL);
props_.viscosity(props_.numCells(), pressure, temperature, NULL, &allcells_[0], &visc_[0], NULL);
props_.matrix(props_.numCells(), pressure, temperature, NULL, &allcells_[0], &A_[0], NULL);
// Check immiscibility requirement (only done for first cell).
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.
/// \param[in] darcyflux Array of signed face fluxes.
/// \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] porevolume0 Array of pore volumes at start 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.
void solve(const double* darcyflux,
const double* pressure,
const double* temperature,
const double* porevolume0,
const double* porevolume,
const double* source,

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

@ -636,7 +636,7 @@ namespace Opm
double mob[max_np];
props.relperm(1, &s[np*cell], &cell, mob, 0);
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;
for(int i = 0; i < np; ++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.");
}
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>& z = state.surfacevol();
std::fill(injected, injected + np, 0.0);
@ -94,8 +95,8 @@ namespace Opm
const double flux = -transport_src[c]*dt;
const double* sat = &s[np*c];
props.relperm(1, sat, &c, &mob[0], 0);
props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0);
props.matrix(1, &press[c], &z[np*c], &c, &A[0], 0);
props.viscosity(1, &press[c], &temp[c], &z[np*c], &c, &visc[0], 0);
props.matrix(1, &press[c], &temp[c], &z[np*c], &c, &A[0], 0);
double totmob = 0.0;
for (int p = 0; p < np; ++p) {
mob[p] /= visc[p];
@ -121,19 +122,21 @@ namespace Opm
/// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated
/// @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] s saturation values (for all phases)
/// @param[out] totmob total mobilities.
void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells,
const std::vector<double>& press,
const std::vector<double>& temp,
const std::vector<double>& z,
const std::vector<double>& s,
std::vector<double>& totmob)
{
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::vector<int>::size_type nc = cells.size();
@ -193,12 +196,14 @@ namespace Opm
/// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated
/// @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] s saturation values (for all phases)
/// @param[out] pmobc phase mobilities (for all phases).
void computePhaseMobilities(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& T,
const std::vector<double>& z,
const std::vector<double>& s,
std::vector<double>& pmobc)
@ -209,7 +214,7 @@ namespace Opm
assert(int(s.size()) == 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.resize(nc*np, 0.0);
@ -227,19 +232,21 @@ namespace Opm
/// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated
/// @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] s saturation values (for all phases)
/// @param[out] fractional_flow the fractional flow for each phase for each cell.
void computeFractionalFlow(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& T,
const std::vector<double>& z,
const std::vector<double>& s,
std::vector<double>& fractional_flows)
{
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) {
double phase_sum = 0.0;
@ -299,7 +306,7 @@ namespace Opm
//std::vector<double> res_vol(np);
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.
MAT_SIZE_T n = np;
@ -388,7 +395,7 @@ namespace Opm
} else {
assert(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
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.
}
}

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

@ -92,12 +92,14 @@ namespace Opm
/// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated
/// @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] s saturation values (for all phases)
/// @param[out] pmobc phase mobilities (for all phases).
void computePhaseMobilities(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& T,
const std::vector<double>& z,
const std::vector<double>& s,
std::vector<double>& pmobc);
@ -107,12 +109,14 @@ namespace Opm
/// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated
/// @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] s saturation values (for all phases)
/// @param[out] fractional_flow the fractional flow for each phase for each cell.
void computeFractionalFlow(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& T,
const std::vector<double>& z,
const std::vector<double>& s,
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_;
}
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<double> mu(n);
@ -116,8 +116,8 @@ void testmu(const double reltol, int n, int np, const std::vector<int> &pvtTable
// test mu
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], &z[0], &mu[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], &T[0], &z[0], &mu[0]);
dmudp_diff = (mu_new[1]-mu_new[0])/(p[1]-p[0]);
dmudr_diff = (mu_new[2]-mu_new[0])/(r[2]-r[0]);
dmudp_diff_u = (mu_new[4]-mu_new[3])/(p[4]-p[3]);
@ -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)
{
// test b
@ -155,8 +155,8 @@ void testb(const double reltol, int n, int np, const std::vector<int> &pvtTableI
double dbdr_diff_u;
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]->dBdp(n, &pvtTableIdx[0], &p[0], &z[0], &B[0], &dBdp[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], &T[0], &z[0], &B[0], &dBdp[0]);
dbdp_diff = (b[1]-b[0])/(p[1]-p[0]);
dbdr_diff = (b[2]-b[0])/(r[2]-r[0]);
dbdp_diff_u = (b[4]-b[3])/(p[4]-p[3]);
@ -253,6 +253,7 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
const double reltolpermu = 1e-1;
std::vector<double> p(n);
std::vector<double> T(n, 273.15 + 20);
std::vector<double> r(n);
std::vector<PhasePresence> condition(n);
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_);
@ -332,6 +333,7 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
const double reltolpermu = 1e-1;
std::vector<double> p(n);
std::vector<double> T(n, 273.15+20);
std::vector<double> r(n);
std::vector<PhasePresence> condition(n);
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_);