diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 24ac6438..1ef81767 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -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) { diff --git a/opm/core/props/BlackoilPropertiesBasic.cpp b/opm/core/props/BlackoilPropertiesBasic.cpp index 7ca21ca2..ac690ccd 100644 --- a/opm/core/props/BlackoilPropertiesBasic.cpp +++ b/opm/core/props/BlackoilPropertiesBasic.cpp @@ -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) { diff --git a/opm/core/props/BlackoilPropertiesBasic.hpp b/opm/core/props/BlackoilPropertiesBasic.hpp index 5c0472b9..fcbd19f8 100644 --- a/opm/core/props/BlackoilPropertiesBasic.hpp +++ b/opm/core/props/BlackoilPropertiesBasic.hpp @@ -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, diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index 16cec166..d314849f 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -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(); diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index 3c71bda3..cbcd35a8 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -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, diff --git a/opm/core/props/BlackoilPropertiesInterface.hpp b/opm/core/props/BlackoilPropertiesInterface.hpp index 8a25e344..9ce61ac3 100644 --- a/opm/core/props/BlackoilPropertiesInterface.hpp +++ b/opm/core/props/BlackoilPropertiesInterface.hpp @@ -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, diff --git a/opm/core/props/IncompPropertiesBasic.cpp b/opm/core/props/IncompPropertiesBasic.cpp index 6d7b1093..18a03815 100644 --- a/opm/core/props/IncompPropertiesBasic.cpp +++ b/opm/core/props/IncompPropertiesBasic.cpp @@ -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() diff --git a/opm/core/props/pvt/BlackoilPvtProperties.cpp b/opm/core/props/pvt/BlackoilPvtProperties.cpp index 8e95093e..0cdc456c 100644 --- a/opm/core/props/pvt/BlackoilPvtProperties.cpp +++ b/opm/core/props/pvt/BlackoilPvtProperties.cpp @@ -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]; diff --git a/opm/core/props/pvt/BlackoilPvtProperties.hpp b/opm/core/props/pvt/BlackoilPvtProperties.hpp index 9651dce4..45caa12f 100644 --- a/opm/core/props/pvt/BlackoilPvtProperties.hpp +++ b/opm/core/props/pvt/BlackoilPvtProperties.hpp @@ -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; diff --git a/opm/core/props/pvt/PvtConstCompr.hpp b/opm/core/props/pvt/PvtConstCompr.hpp index ee171ba4..46e68f8f 100644 --- a/opm/core/props/pvt/PvtConstCompr.hpp +++ b/opm/core/props/pvt/PvtConstCompr.hpp @@ -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, diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp index db9cdb39..67ed9335 100644 --- a/opm/core/props/pvt/PvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -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, diff --git a/opm/core/props/pvt/PvtDead.hpp b/opm/core/props/pvt/PvtDead.hpp index ccc9acac..b6f04c2c 100644 --- a/opm/core/props/pvt/PvtDead.hpp +++ b/opm/core/props/pvt/PvtDead.hpp @@ -50,64 +50,71 @@ namespace Opm void initFromGas(const std::vector& 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, diff --git a/opm/core/props/pvt/PvtDeadSpline.cpp b/opm/core/props/pvt/PvtDeadSpline.cpp index 180c0e5a..5d0f635a 100644 --- a/opm/core/props/pvt/PvtDeadSpline.cpp +++ b/opm/core/props/pvt/PvtDeadSpline.cpp @@ -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, diff --git a/opm/core/props/pvt/PvtDeadSpline.hpp b/opm/core/props/pvt/PvtDeadSpline.hpp index 83b14c17..3f60f9f3 100644 --- a/opm/core/props/pvt/PvtDeadSpline.hpp +++ b/opm/core/props/pvt/PvtDeadSpline.hpp @@ -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, diff --git a/opm/core/props/pvt/PvtInterface.hpp b/opm/core/props/pvt/PvtInterface.hpp index 85b5b74f..3f25f21b 100644 --- a/opm/core/props/pvt/PvtInterface.hpp +++ b/opm/core/props/pvt/PvtInterface.hpp @@ -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, diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp index 72520613..ce189c8b 100644 --- a/opm/core/props/pvt/PvtLiveGas.cpp +++ b/opm/core/props/pvt/PvtLiveGas.cpp @@ -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, diff --git a/opm/core/props/pvt/PvtLiveGas.hpp b/opm/core/props/pvt/PvtLiveGas.hpp index 945502fe..a54a793f 100644 --- a/opm/core/props/pvt/PvtLiveGas.hpp +++ b/opm/core/props/pvt/PvtLiveGas.hpp @@ -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& 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, diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp index 012a1fdb..bf78e48f 100644 --- a/opm/core/props/pvt/PvtLiveOil.cpp +++ b/opm/core/props/pvt/PvtLiveOil.cpp @@ -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, diff --git a/opm/core/props/pvt/PvtLiveOil.hpp b/opm/core/props/pvt/PvtLiveOil.hpp index 2b3a69dc..e94dd7c1 100644 --- a/opm/core/props/pvt/PvtLiveOil.hpp +++ b/opm/core/props/pvt/PvtLiveOil.hpp @@ -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& 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, diff --git a/opm/core/props/pvt/PvtPropertiesBasic.cpp b/opm/core/props/pvt/PvtPropertiesBasic.cpp index 29790522..d136f12e 100644 --- a/opm/core/props/pvt/PvtPropertiesBasic.cpp +++ b/opm/core/props/pvt/PvtPropertiesBasic.cpp @@ -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 diff --git a/opm/core/props/pvt/PvtPropertiesBasic.hpp b/opm/core/props/pvt/PvtPropertiesBasic.hpp index 46add622..ca444f9f 100644 --- a/opm/core/props/pvt/PvtPropertiesBasic.hpp +++ b/opm/core/props/pvt/PvtPropertiesBasic.hpp @@ -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; diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index f7c76435..16095f44 100644 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ b/opm/core/simulator/EquilibrationHelpers.hpp @@ -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 operator()(const double p, + const double T, const std::vector& z) const { const int np = props_.numPhases(); @@ -136,7 +139,7 @@ namespace Opm assert (z.size() == std::vector::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 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. diff --git a/opm/core/simulator/SimulatorCompressibleTwophase.cpp b/opm/core/simulator/SimulatorCompressibleTwophase.cpp index 7284d8d4..730306e2 100644 --- a/opm/core/simulator/SimulatorCompressibleTwophase.cpp +++ b/opm/core/simulator/SimulatorCompressibleTwophase.cpp @@ -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 }; diff --git a/opm/core/simulator/SimulatorState.cpp b/opm/core/simulator/SimulatorState.cpp index b02df8bd..6b2eac57 100644 --- a/opm/core/simulator/SimulatorState.cpp +++ b/opm/core/simulator/SimulatorState.cpp @@ -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); diff --git a/opm/core/simulator/SimulatorState.hpp b/opm/core/simulator/SimulatorState.hpp index 5f4fbf08..71a0a605 100644 --- a/opm/core/simulator/SimulatorState.hpp +++ b/opm/core/simulator/SimulatorState.hpp @@ -38,11 +38,13 @@ namespace Opm int numPhases() const { return num_phases_; } std::vector& pressure () { return press_ ; } + std::vector& temperature () { return temp_ ; } std::vector& facepressure() { return fpress_; } std::vector& faceflux () { return flux_ ; } std::vector& saturation () { return sat_ ; } const std::vector& pressure () const { return press_ ; } + const std::vector& temperature () const { return temp_ ; } const std::vector& facepressure() const { return fpress_; } const std::vector& faceflux () const { return flux_ ; } const std::vector& saturation () const { return sat_ ; } @@ -56,6 +58,7 @@ namespace Opm private: int num_phases_; std::vector press_ ; + std::vector temp_ ; std::vector fpress_; std::vector flux_ ; std::vector sat_ ; diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 090676ee..02fa09f0 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -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& bhp() { return bhp_; } const std::vector& bhp() const { return bhp_; } + /// One temperature per well. + std::vector& temperature() { return temperature_; } + const std::vector& temperature() const { return temperature_; } + /// One rate per well and phase. std::vector& wellRates() { return wellrates_; } const std::vector& wellRates() const { return wellrates_; } @@ -124,6 +129,7 @@ namespace Opm private: std::vector bhp_; + std::vector temperature_; std::vector wellrates_; std::vector perfrates_; std::vector perfpress_; diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 07c34129..8e86c4ef 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -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 computeRs(const UnstructuredGrid& grid, const CellRangeType& cells, const std::vector oil_pressure, + const std::vector& temperature, const Miscibility::RsFunction& rs_func, const std::vector 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(props, cell, p_contact)); + const double T_contact = 273.15 + 20; // standard temperature for now + rs_func_.push_back(std::make_shared(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(props, cell, p_contact)); + const double T_contact = 273.15 + 20; // standard temperature for now + rv_func_.push_back(std::make_shared(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& 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_); } diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index 9656d13c..b31297b8 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -106,11 +106,13 @@ namespace Opm template 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 svol_; const int ix_; @@ -134,7 +137,7 @@ namespace Opm double density(const double press) const { - const std::vector& rho = rho_(press, svol_); + const std::vector& rho = rho_(press, temp_, svol_); return rho[ix_]; } @@ -143,13 +146,15 @@ namespace Opm template 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 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& rho = rho_(press, svol_); + const std::vector& rho = rho_(press, temp_, svol_); return rho[oix_]; } }; @@ -190,13 +196,15 @@ namespace Opm template 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 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& rho = rho_(press, svol_); + const std::vector& 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 + std::vector + temperature(const UnstructuredGrid& G, + const Region& reg, + const CellRange& cells) + { + // use the standard temperature for everything for now + return std::vector(cells.size(), 273.15 + 20.0); + } template std::vector< std::vector > @@ -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 computeRs(const UnstructuredGrid& grid, const CellRangeType& cells, const std::vector oil_pressure, + const std::vector& temperature, const Miscibility::RsFunction& rs_func, const std::vector 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; } diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index 70a6248d..1879a215 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -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 diff --git a/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp index 82f8ea80..278e500c 100644 --- a/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp +++ b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp @@ -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) { diff --git a/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp index c1c9f3fb..4092eb11 100644 --- a/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp +++ b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp @@ -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, diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 4b96ead7..b7f409b8 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -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]; diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 8a0ebac8..f650b945 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -74,6 +74,7 @@ namespace Opm OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells."); } const std::vector& press = state.pressure(); + const std::vector& temp = state.temperature(); const std::vector& s = state.saturation(); const std::vector& 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& cells, const std::vector& press, + const std::vector& temp, const std::vector& z, const std::vector& s, std::vector& totmob) { std::vector 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::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& cells, const std::vector& p, + const std::vector& T, const std::vector& z, const std::vector& s, std::vector& pmobc) @@ -209,7 +214,7 @@ namespace Opm assert(int(s.size()) == nc * np); std::vector 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& cells, const std::vector& p, + const std::vector& T, const std::vector& z, const std::vector& s, std::vector& 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::size_type i = 0; i < cells.size(); ++i) { double phase_sum = 0.0; @@ -299,7 +306,7 @@ namespace Opm //std::vector res_vol(np); const std::vector& 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. } } diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp index 40d59e75..786a5c30 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.hpp +++ b/opm/core/utility/miscUtilitiesBlackoil.hpp @@ -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& cells, const std::vector& p, + const std::vector& T, const std::vector& z, const std::vector& s, std::vector& 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& cells, const std::vector& p, + const std::vector& T, const std::vector& z, const std::vector& s, std::vector& fractional_flows); diff --git a/tests/test_blackoilfluid.cpp b/tests/test_blackoilfluid.cpp index add7d495..aed8b45d 100644 --- a/tests/test_blackoilfluid.cpp +++ b/tests/test_blackoilfluid.cpp @@ -102,7 +102,7 @@ std::vector > getProps(Opm::DeckConstPtr deck, Opm return props_; } -void testmu(const double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector r,std::vector z, +void testmu(const double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector T, std::vector r,std::vector z, std::vector > props_, std::vector condition) { std::vector mu(n); @@ -116,8 +116,8 @@ void testmu(const double reltol, int n, int np, const std::vector &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 &pvtTable } } -void testb(const double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector r,std::vector z, +void testb(const double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector T, std::vector r,std::vector z, std::vector > props_, std::vector condition) { // test b @@ -155,8 +155,8 @@ void testb(const double reltol, int n, int np, const std::vector &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 p(n); + std::vector T(n, 273.15 + 20); std::vector r(n); std::vector condition(n); std::vector 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 p(n); + std::vector T(n, 273.15+20); std::vector r(n); std::vector condition(n); std::vector 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_);