From e530ce03d8b878f7834976408f7cf8fee079f2e2 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 20 Nov 2014 12:15:01 +0100 Subject: [PATCH] PVT properties: allow them to be temperature dependent 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... --- opm/core/pressure/CompressibleTpfa.cpp | 12 ++-- opm/core/props/BlackoilPropertiesBasic.cpp | 10 ++- opm/core/props/BlackoilPropertiesBasic.hpp | 4 ++ opm/core/props/BlackoilPropertiesFromDeck.cpp | 10 ++- opm/core/props/BlackoilPropertiesFromDeck.hpp | 4 ++ .../props/BlackoilPropertiesInterface.hpp | 4 ++ opm/core/props/IncompPropertiesBasic.cpp | 4 +- opm/core/props/pvt/BlackoilPvtProperties.cpp | 9 ++- opm/core/props/pvt/BlackoilPvtProperties.hpp | 9 ++- opm/core/props/pvt/PvtConstCompr.hpp | 7 ++ opm/core/props/pvt/PvtDead.cpp | 23 +++--- opm/core/props/pvt/PvtDead.hpp | 21 ++++-- opm/core/props/pvt/PvtDeadSpline.cpp | 15 ++-- opm/core/props/pvt/PvtDeadSpline.hpp | 21 ++++-- opm/core/props/pvt/PvtInterface.hpp | 31 ++++---- opm/core/props/pvt/PvtLiveGas.cpp | 31 ++++---- opm/core/props/pvt/PvtLiveGas.hpp | 27 ++++--- opm/core/props/pvt/PvtLiveOil.cpp | 37 ++++++---- opm/core/props/pvt/PvtLiveOil.hpp | 25 ++++--- opm/core/props/pvt/PvtPropertiesBasic.cpp | 3 + opm/core/props/pvt/PvtPropertiesBasic.hpp | 11 +-- opm/core/simulator/EquilibrationHelpers.hpp | 71 +++++++++++++------ .../SimulatorCompressibleTwophase.cpp | 4 +- opm/core/simulator/SimulatorState.cpp | 2 + opm/core/simulator/SimulatorState.hpp | 3 + opm/core/simulator/WellState.hpp | 6 ++ opm/core/simulator/initStateEquil.hpp | 13 ++-- opm/core/simulator/initStateEquil_impl.hpp | 58 ++++++++++----- opm/core/simulator/initState_impl.hpp | 19 ++--- ...sportSolverCompressibleTwophaseReorder.cpp | 5 +- ...sportSolverCompressibleTwophaseReorder.hpp | 2 + opm/core/utility/miscUtilities.cpp | 2 +- opm/core/utility/miscUtilitiesBlackoil.cpp | 21 ++++-- opm/core/utility/miscUtilitiesBlackoil.hpp | 4 ++ tests/test_blackoilfluid.cpp | 22 +++--- 35 files changed, 370 insertions(+), 180 deletions(-) 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_);