From 0aa7620a0e34a8df795138d04c3413ea42e4f777 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 17 Mar 2015 12:40:05 +0100 Subject: [PATCH 1/4] PVT classes: do not implement temperature dependence directly anymore instead, we will be going with wrapper classes around the PvtInterface in the next commits. this considerably reduces the amount of copy-and-paste required for temperature support. --- opm/core/props/pvt/PvtConstCompr.hpp | 77 +--------------------------- opm/core/props/pvt/PvtDead.cpp | 46 +---------------- opm/core/props/pvt/PvtDead.hpp | 18 +------ opm/core/props/pvt/PvtDeadSpline.cpp | 43 ++-------------- opm/core/props/pvt/PvtDeadSpline.hpp | 12 ----- opm/core/props/pvt/PvtLiveGas.cpp | 5 -- opm/core/props/pvt/PvtLiveOil.cpp | 51 ++---------------- opm/core/props/pvt/PvtLiveOil.hpp | 11 ---- 8 files changed, 11 insertions(+), 252 deletions(-) diff --git a/opm/core/props/pvt/PvtConstCompr.hpp b/opm/core/props/pvt/PvtConstCompr.hpp index 1e76f97c..46e68f8f 100644 --- a/opm/core/props/pvt/PvtConstCompr.hpp +++ b/opm/core/props/pvt/PvtConstCompr.hpp @@ -23,8 +23,6 @@ #include #include -#include - #include #include @@ -47,8 +45,6 @@ namespace Opm { public: PvtConstCompr() - : oilvisctTables_(0), - watvisctTables_(0) {} void initFromWater(Opm::DeckKeywordConstPtr pvtwKeyword) @@ -102,9 +98,7 @@ namespace Opm ref_B_(1, 1.0), comp_(1, 0.0), viscosity_(1, visc), - visc_comp_(1, 0.0), - oilvisctTables_(0), - watvisctTables_(0) + visc_comp_(1, 0.0) { } @@ -125,13 +119,6 @@ namespace Opm int tableIdx = getTableIndex_(pvtRegionIdx, i); double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); output_mu[i] = viscosity_[tableIdx]/(1.0 + x + 0.5*x*x); - - if (oilvisctTables_ != 0 || watvisctTables_ != 0) { - // TODO: temperature dependence - OPM_THROW(std::logic_error, - "temperature dependent viscosity as a function of z " - "is not yet implemented!"); - } } } @@ -159,7 +146,7 @@ namespace Opm virtual void mu(const int n, const int* pvtRegionIdx, const double* p, - const double* T, + const double* /*T*/, const double* /*r*/, const PhasePresence* /*cond*/, double* output_mu, @@ -174,46 +161,6 @@ namespace Opm double d = (1.0 + x + 0.5*x*x); output_mu[i] = viscosity_[tableIdx]/d; output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx]; - - if (oilvisctTables_ != 0) { - // this object handles _either_ oil _or_ water - assert(watvisctTables_ == 0); - - // temperature dependence of the oil phase - DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx); - double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); - double muRef = -visc_comp_[tableIdx]*(pRef - ref_press_[tableIdx]); - - double muOilvisct = (*oilvisctTables_)[tableIdx].evaluate("Viscosity", T[i]); - double alpha = muOilvisct/muRef; - - output_mu[i] *= alpha; - output_dmudp[i] *= alpha; - output_dmudr[i] *= alpha; - - // TODO (?): derivative of oil viscosity w.r.t. temperature. - // probably requires a healthy portion of if-spaghetti - } - - if (watvisctTables_ != 0) { - // this object handles _either_ oil _or_ water - assert(oilvisctTables_ == 0); - - // temperature dependence of the water phase - DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx); - double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); - double muRef = -visc_comp_[tableIdx]*(pRef - ref_press_[tableIdx]); - - double muWatvisct = (*watvisctTables_)[tableIdx].evaluate("Viscosity", T[i]); - double alpha = muWatvisct/muRef; - - output_mu[i] *= alpha; - output_dmudp[i] *= alpha; - output_dmudr[i] *= alpha; - - // TODO (?): derivative of oil viscosity w.r.t. temperature. - // probably requires a healthy portion of if-spaghetti - } } std::fill(output_dmudr, output_dmudr + n, 0.0); } @@ -340,22 +287,6 @@ namespace Opm std::fill(output_dRdp, output_dRdp + n, 0.0); } - /// set the tables which specify the temperature dependence of the oil viscosity - void setOilvisctTables(const std::vector& oilvisctTables, - DeckKeywordConstPtr viscrefKeyword) - { - oilvisctTables_ = &oilvisctTables; - viscrefKeyword_ = viscrefKeyword; - } - - /// set the tables which specify the temperature dependence of the water viscosity - void setWatvisctTables(const std::vector& watvisctTables, - DeckKeywordConstPtr viscrefKeyword) - { - watvisctTables_ = &watvisctTables; - viscrefKeyword_ = viscrefKeyword; - } - private: int getTableIndex_(const int* pvtTableIdx, int cellIdx) const { @@ -371,10 +302,6 @@ namespace Opm std::vector comp_; std::vector viscosity_; std::vector visc_comp_; - - const std::vector* oilvisctTables_; - const std::vector* watvisctTables_; - DeckKeywordConstPtr viscrefKeyword_; }; } diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp index 3d711ca6..67ed9335 100644 --- a/opm/core/props/pvt/PvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -123,20 +123,13 @@ namespace Opm double tempInvB = b_[regionIdx](p[i]); double tempInvBmu = inverseBmu_[regionIdx](p[i]); output_mu[i] = tempInvB / tempInvBmu; - - if (oilvisctTables_ != 0) { - // TODO: temperature dependence - OPM_THROW(std::logic_error, - "temperature dependent viscosity as a function of z " - "is not yet implemented!"); - } } } void PvtDead::mu(const int n, const int* pvtTableIdx, const double* p, - const double* T, + const double* /*T*/, const double* /*r*/, double* output_mu, double* output_dmudp, @@ -150,24 +143,6 @@ namespace Opm output_mu[i] = tempInvB / tempInvBmu; output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i]) - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu); - - if (oilvisctTables_ != 0) { - // temperature dependence of the oil phase - DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx); - double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); - double muRef = b_[regionIdx](pRef)/inverseBmu_[regionIdx](pRef); - - double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]); - double alpha = muOilvisct/muRef; - - output_mu[i] *= alpha; - output_dmudp[i] *= alpha; - output_dmudr[i] *= alpha; - - // TODO (?): derivative of oil viscosity w.r.t. temperature. - // probably requires a healthy portion of if-spaghetti - } - } std::fill(output_dmudr, output_dmudr + n, 0.0); @@ -176,7 +151,7 @@ namespace Opm void PvtDead::mu(const int n, const int* pvtTableIdx, const double* p, - const double* T, + const double* /*T*/, const double* /*r*/, const PhasePresence* /*cond*/, double* output_mu, @@ -192,23 +167,6 @@ namespace Opm output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i]) - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu); - - if (oilvisctTables_ != 0) { - // temperature dependence of the oil phase - DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx); - double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); - double muRef = b_[regionIdx](pRef)/inverseBmu_[regionIdx](pRef); - - double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]); - double alpha = muOilvisct/muRef; - - output_mu[i] *= alpha; - output_dmudp[i] *= alpha; - output_dmudr[i] *= alpha; - - // TODO (?): derivative of oil viscosity w.r.t. temperature. - // probably requires a healthy portion of if-spaghetti - } } std::fill(output_dmudr, output_dmudr + n, 0.0); diff --git a/opm/core/props/pvt/PvtDead.hpp b/opm/core/props/pvt/PvtDead.hpp index 0bb20ea6..b6f04c2c 100644 --- a/opm/core/props/pvt/PvtDead.hpp +++ b/opm/core/props/pvt/PvtDead.hpp @@ -44,11 +44,7 @@ namespace Opm class PvtDead : public PvtInterface { public: - PvtDead() - { - // by default, specify no temperature dependence of the PVT properties - oilvisctTables_ = 0; - } + PvtDead() {}; void initFromOil(const std::vector& pvdoTables); void initFromGas(const std::vector& pvdgTables); @@ -154,15 +150,6 @@ namespace Opm const double* z, double* output_R, double* output_dRdp) const; - - /// set the tables which specify the temperature dependence of the oil viscosity - void setOilvisctTables(const std::vector& oilvisctTables, - DeckKeywordConstPtr viscrefKeyword) - { - oilvisctTables_ = &oilvisctTables; - viscrefKeyword_ = viscrefKeyword; - } - private: int getTableIndex_(const int* pvtTableIdx, int cellIdx) const { @@ -176,9 +163,6 @@ namespace Opm std::vector > b_; std::vector > viscosity_; std::vector > inverseBmu_; - - const std::vector* oilvisctTables_; - DeckKeywordConstPtr viscrefKeyword_; }; } diff --git a/opm/core/props/pvt/PvtDeadSpline.cpp b/opm/core/props/pvt/PvtDeadSpline.cpp index 1e0b37e9..293caf62 100644 --- a/opm/core/props/pvt/PvtDeadSpline.cpp +++ b/opm/core/props/pvt/PvtDeadSpline.cpp @@ -36,10 +36,7 @@ namespace Opm //------------------------------------------------------------------------- PvtDeadSpline::PvtDeadSpline() - { - // by default, specify no temperature dependence of the PVT properties - oilvisctTables_ = 0; - } + {} void PvtDeadSpline::initFromOil(const std::vector& pvdoTables, int numSamples) @@ -123,7 +120,7 @@ namespace Opm void PvtDeadSpline::mu(const int n, const int* pvtTableIdx, const double* p, - const double* T, + const double* /*T*/, const double* /*r*/, double* output_mu, double* output_dmudp, @@ -134,23 +131,6 @@ namespace Opm int regionIdx = getTableIndex_(pvtTableIdx, i); output_mu[i] = viscosity_[regionIdx](p[i]); output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); - - if (oilvisctTables_ != 0) { - // temperature dependence of the oil phase - DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx); - double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); - double muRef = viscosity_[regionIdx](pRef); - - double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]); - double alpha = muOilvisct/muRef; - - output_mu[i] *= alpha; - output_dmudp[i] *= alpha; - output_dmudr[i] *= alpha; - - // TODO (?): derivative of oil viscosity w.r.t. temperature. - // probably requires a healthy portion of if-spaghetti - } } std::fill(output_dmudr, output_dmudr + n, 0.0); } @@ -158,7 +138,7 @@ namespace Opm void PvtDeadSpline::mu(const int n, const int* pvtTableIdx, const double* p, - const double* T, + const double* /*T*/, const double* /*r*/, const PhasePresence* /*cond*/, double* output_mu, @@ -171,23 +151,6 @@ namespace Opm int regionIdx = getTableIndex_(pvtTableIdx, i); output_mu[i] = viscosity_[regionIdx](p[i]); output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); - - if (oilvisctTables_ != 0) { - // temperature dependence of the oil phase - DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx); - double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); - double muRef = viscosity_[regionIdx](pRef); - - double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]); - double alpha = muOilvisct/muRef; - - output_mu[i] *= alpha; - output_dmudp[i] *= alpha; - output_dmudr[i] *= alpha; - - // TODO (?): derivative of oil viscosity w.r.t. temperature. - // probably requires a healthy portion of if-spaghetti - } } std::fill(output_dmudr, output_dmudr + n, 0.0); } diff --git a/opm/core/props/pvt/PvtDeadSpline.hpp b/opm/core/props/pvt/PvtDeadSpline.hpp index 366e0caf..3f60f9f3 100644 --- a/opm/core/props/pvt/PvtDeadSpline.hpp +++ b/opm/core/props/pvt/PvtDeadSpline.hpp @@ -148,15 +148,6 @@ namespace Opm const double* z, double* output_R, double* output_dRdp) const; - - /// set the tables which specify the temperature dependence of the oil viscosity - void setOilvisctTables(const std::vector& oilvisctTables, - DeckKeywordConstPtr viscrefKeyword) - { - oilvisctTables_ = &oilvisctTables; - viscrefKeyword_ = viscrefKeyword; - } - private: int getTableIndex_(const int* pvtTableIdx, int cellIdx) const { @@ -169,9 +160,6 @@ namespace Opm // table per PVT region. std::vector > b_; std::vector > viscosity_; - - const std::vector* oilvisctTables_; - DeckKeywordConstPtr viscrefKeyword_; }; } diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp index cf60b836..ce189c8b 100644 --- a/opm/core/props/pvt/PvtLiveGas.cpp +++ b/opm/core/props/pvt/PvtLiveGas.cpp @@ -114,9 +114,6 @@ namespace Opm double inverseBMu = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 3, false); output_mu[i] = inverseB / inverseBMu; - - // temperature dependence: Since E100 does not implement temperature - // dependence of gas viscosity, we skip it here as well... } } @@ -165,8 +162,6 @@ namespace Opm output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr) / (inverseBMu * inverseBMu); - // temperature dependence: Since E100 does not implement temperature - // dependence of gas viscosity, we skip it here as well... } } diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp index 9d28c0fa..bf78e48f 100644 --- a/opm/core/props/pvt/PvtLiveOil.cpp +++ b/opm/core/props/pvt/PvtLiveOil.cpp @@ -36,9 +36,6 @@ namespace Opm //------------------------------------------------------------------------- PvtLiveOil::PvtLiveOil(const std::vector& pvtoTables) { - // by default, specify no temperature dependence of the PVT properties - oilvisctTables_ = 0; - int numTables = pvtoTables.size(); saturated_oil_table_.resize(numTables); undersat_oil_tables_.resize(numTables); @@ -152,20 +149,13 @@ namespace Opm output_mu[i] = inverseB / inverseBMu; } - - if (oilvisctTables_ != 0) { - // TODO: temperature dependence - OPM_THROW(std::logic_error, - "temperature dependent viscosity as a function of z " - "is not yet implemented!"); - } } /// 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* T, + const double* /*T*/, const double* r, double* output_mu, double* output_dmudp, @@ -193,32 +183,14 @@ namespace Opm output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr) / (inverseBMu * inverseBMu); - if (oilvisctTables_ != 0) { - // temperature dependence - DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx); - double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); - double rsRef = viscrefRecord->getItem("REFERENCE_RS")->getSIDouble(0); - double muRef = miscible_oil(pRef, rsRef, tableIdx, 3, 0)/miscible_oil(pRef, rsRef, tableIdx, 1, 0); - - double muOilvisct = (*oilvisctTables_)[tableIdx].evaluate("Viscosity", T[i]); - double alpha = muOilvisct/muRef; - - output_mu[i] *= alpha; - output_dmudp[i] *= alpha; - output_dmudr[i] *= alpha; - - // TODO (?): derivative of oil viscosity w.r.t. temperature. - // probably requires a healthy portion of if-spaghetti - } } - } /// 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* T, + const double* /*T*/, const double* r, const PhasePresence* cond, double* output_mu, @@ -248,24 +220,6 @@ namespace Opm output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr) / (inverseBMu * inverseBMu); - - if (oilvisctTables_ != 0) { - // temperature dependence - DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx); - double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); - double rsRef = viscrefRecord->getItem("REFERENCE_RS")->getSIDouble(0); - double muRef = miscible_oil(pRef, rsRef, tableIdx, 3, 0)/miscible_oil(pRef, rsRef, tableIdx, 1, 0); - - double muOilvisct = (*oilvisctTables_)[tableIdx].evaluate("Viscosity", T[i]); - double alpha = muOilvisct/muRef; - - output_mu[i] *= alpha; - output_dmudp[i] *= alpha; - output_dmudr[i] *= alpha; - - // TODO (?): derivative of oil viscosity w.r.t. temperature. - // probably requires a healthy portion of if-spaghetti - } } } @@ -691,4 +645,5 @@ namespace Opm } } } + } // namespace Opm diff --git a/opm/core/props/pvt/PvtLiveOil.hpp b/opm/core/props/pvt/PvtLiveOil.hpp index 34191875..e94dd7c1 100644 --- a/opm/core/props/pvt/PvtLiveOil.hpp +++ b/opm/core/props/pvt/PvtLiveOil.hpp @@ -142,14 +142,6 @@ namespace Opm double* output_R, double* output_dRdp) const; - /// set the tables which specify the temperature dependence of the oil viscosity - void setOilvisctTables(const std::vector& oilvisctTables, - DeckKeywordConstPtr viscrefKeyword) - { - oilvisctTables_ = &oilvisctTables; - viscrefKeyword_ = viscrefKeyword; - } - private: int getTableIndex_(const int* pvtTableIdx, int cellIdx) const { @@ -187,9 +179,6 @@ namespace Opm // store one table per PVT region. std::vector > > saturated_oil_table_; std::vector > > > undersat_oil_tables_; - - const std::vector* oilvisctTables_; - DeckKeywordConstPtr viscrefKeyword_; }; } From 2a52c2926d886ab59a63348853958933aa1e880a Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 17 Mar 2015 12:40:07 +0100 Subject: [PATCH 2/4] add a wrapper for thermal water PVT objects --- CMakeLists_files.cmake | 1 + opm/core/props/pvt/ThermalWaterPvtWrapper.hpp | 390 ++++++++++++++++++ 2 files changed, 391 insertions(+) create mode 100644 opm/core/props/pvt/ThermalWaterPvtWrapper.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index af66a669..5a4d37a5 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -356,6 +356,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/props/pvt/PvtInterface.hpp opm/core/props/pvt/PvtLiveGas.hpp opm/core/props/pvt/PvtLiveOil.hpp + opm/core/props/pvt/ThermalWaterPvtWrapper.hpp opm/core/props/rock/RockBasic.hpp opm/core/props/rock/RockCompressibility.hpp opm/core/props/rock/RockFromDeck.hpp diff --git a/opm/core/props/pvt/ThermalWaterPvtWrapper.hpp b/opm/core/props/pvt/ThermalWaterPvtWrapper.hpp new file mode 100644 index 00000000..baf46c4d --- /dev/null +++ b/opm/core/props/pvt/ThermalWaterPvtWrapper.hpp @@ -0,0 +1,390 @@ +/* + Copyright 2015 Andreas Lauser + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_THERMAL_WATER_PVT_WRAPPER_HPP +#define OPM_THERMAL_WATER_PVT_WRAPPER_HPP + +#include +#include + +#include + +#include + +namespace Opm +{ + /// Class which wraps another (i.e., isothermal) PVT object into one which adds + /// temperature dependence of water + class ThermalWaterPvtWrapper : public PvtInterface + { + public: + ThermalWaterPvtWrapper() + {} + + + /// set the tables which specify the temperature dependence of the water viscosity + void initFromDeck(std::shared_ptr isothermalPvt, + Opm::DeckConstPtr deck, + Opm::EclipseStateConstPtr eclipseState) + { + isothermalPvt_ = isothermalPvt; + watvisctTables_ = 0; + + // stuff which we need to get from the PVTW keyword + Opm::DeckKeywordConstPtr pvtwKeyword = deck->getKeyword("PVTW"); + int numRegions = pvtwKeyword->size(); + pvtwRefPress_.resize(numRegions); + pvtwRefB_.resize(numRegions); + pvtwCompressibility_.resize(numRegions); + pvtwViscosity_.resize(numRegions); + pvtwViscosibility_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) { + Opm::DeckRecordConstPtr pvtwRecord = pvtwKeyword->getRecord(regionIdx); + pvtwRefPress_[regionIdx] = pvtwRecord->getItem("P_REF")->getSIDouble(0); + pvtwRefB_[regionIdx] = pvtwRecord->getItem("WATER_VOL_FACTOR")->getSIDouble(0); + pvtwViscosity_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSITY")->getSIDouble(0); + pvtwViscosibility_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSIBILITY")->getSIDouble(0); + } + + // quantities required for the temperature dependence of the viscosity + // (basically we expect well-behaved VISCREF and WATVISCT keywords.) + if (deck->hasKeyword("VISCREF")) { + watvisctTables_ = &eclipseState->getWatvisctTables(); + Opm::DeckKeywordConstPtr viscrefKeyword = deck->getKeyword("VISCREF"); + + assert(watvisctTables_->size() == numRegions); + assert(viscrefKeyword->size() == numRegions); + + viscrefPress_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) { + Opm::DeckRecordConstPtr viscrefRecord = viscrefKeyword->getRecord(regionIdx); + + viscrefPress_[regionIdx] = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); + } + } + + // quantities required for the temperature dependence of the density + if (deck->hasKeyword("WATDENT")) { + DeckKeywordConstPtr watdentKeyword = deck->getKeyword("WATDENT"); + + assert(watdentKeyword->size() == numRegions); + + watdentRefTemp_.resize(numRegions); + watdentCT1_.resize(numRegions); + watdentCT2_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + Opm::DeckRecordConstPtr watdentRecord = watdentKeyword->getRecord(regionIdx); + + watdentRefTemp_[regionIdx] = watdentRecord->getItem("REFERENCE_TEMPERATURE")->getSIDouble(0); + watdentCT1_[regionIdx] = watdentRecord->getItem("EXPANSION_COEFF_LINEAR")->getSIDouble(0); + watdentCT2_[regionIdx] = watdentRecord->getItem("EXPANSION_COEFF_QUADRATIC")->getSIDouble(0); + } + } + } + + virtual void mu(const int n, + const int* pvtRegionIdx, + const double* p, + const double* T, + const double* z, + double* output_mu) const + { + if (watvisctTables_) + // TODO: temperature dependence for viscosity depending on z + OPM_THROW(std::runtime_error, + "temperature dependent viscosity as a function of z " + "is not yet implemented!"); + + // compute the isothermal viscosity + isothermalPvt_->mu(n, pvtRegionIdx, p, T, z, output_mu); + } + + 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 + { + // compute the isothermal viscosity and its derivatives + isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, output_mu, output_dmudp, output_dmudr); + + if (!watvisctTables_) + // isothermal case + return; + + // temperature dependence + for (int i = 0; i < n; ++i) { + int tableIdx = getTableIndex_(pvtRegionIdx, i); + + // calculate the viscosity of the isothermal keyword for the reference + // pressure given by the VISCREF keyword. + double x = -pvtwViscosibility_[tableIdx]*(viscrefPress_[tableIdx] - pvtwRefPress_[tableIdx]); + double muRef = pvtwViscosity_[tableIdx]/(1.0 + x + 0.5*x*x); + + // compute the viscosity deviation due to temperature + double muWatvisct = (*watvisctTables_)[tableIdx].evaluate("Viscosity", T[i]); + double alpha = muWatvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + // TODO (?): derivative of viscosity w.r.t. temperature. + } + } + + 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 + { + // compute the isothermal viscosity and its derivatives + isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, cond, output_mu, output_dmudp, output_dmudr); + + if (!watvisctTables_) + // isothermal case + return; + + // temperature dependence + for (int i = 0; i < n; ++i) { + int tableIdx = getTableIndex_(pvtRegionIdx, i); + + // calculate the viscosity of the isothermal keyword for the reference + // pressure given by the VISCREF keyword. + double x = -pvtwViscosibility_[tableIdx]*(viscrefPress_[tableIdx] - pvtwRefPress_[tableIdx]); + double muRef = pvtwViscosity_[tableIdx]/(1.0 + x + 0.5*x*x); + + // compute the viscosity deviation due to temperature + double muWatvisct = (*watvisctTables_)[tableIdx].evaluate("Viscosity", T[i]); + double alpha = muWatvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + // TODO (?): derivative of viscosity w.r.t. temperature. + } + } + + virtual void B(const int n, + const int* pvtRegionIdx, + const double* p, + const double* T, + const double* z, + double* output_B) const + { + if (watdentRefTemp_.empty()) { + // isothermal case + isothermalPvt_->B(n, pvtRegionIdx, p, T, z, output_B); + return; + } + + // This changes how the water density depends on pressure compared to what's + // used for the PVTW keyword, but it seems to be what Eclipse does. For + // details, see the documentation for the WATDENT keyword in the Eclipse RM. + for (int i = 0; i < n; ++i) { + int tableIdx = getTableIndex_(pvtRegionIdx, i); + double BwRef = pvtwRefB_[tableIdx]; + double TRef = watdentRefTemp_[tableIdx]; + double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]); + double cT1 = watdentCT1_[tableIdx]; + double cT2 = watdentCT2_[tableIdx]; + double Y = T[i] - TRef; + double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y); + output_B[i] = Bw; + } + } + + 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 + { + if (watdentRefTemp_.empty()) { + // isothermal case + isothermalPvt_->dBdp(n, pvtRegionIdx, p, T, z, output_B, output_dBdp); + return; + } + + // This changes how the water density depends on pressure. This is awkward, + // but it seems to be what Eclipse does. See the documentation for the + // WATDENT keyword in the Eclipse RM + for (int i = 0; i < n; ++i) { + int tableIdx = getTableIndex_(pvtRegionIdx, i); + double BwRef = pvtwRefB_[tableIdx]; + double TRef = watdentRefTemp_[tableIdx]; + double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]); + double cT1 = watdentCT1_[tableIdx]; + double cT2 = watdentCT2_[tableIdx]; + double Y = T[i] - TRef; + double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y); + output_B[i] = Bw; + } + + std::fill(output_dBdp, output_dBdp + n, 0.0); + } + + 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 + { + if (watdentRefTemp_.empty()) { + // isothermal case + isothermalPvt_->b(n, pvtRegionIdx, p, T, r, output_b, output_dbdp, output_dbdr); + return; + } + + // This changes how the water density depends on pressure. This is awkward, + // but it seems to be what Eclipse does. See the documentation for the + // WATDENT keyword in the Eclipse RM + for (int i = 0; i < n; ++i) { + int tableIdx = getTableIndex_(pvtRegionIdx, i); + double BwRef = pvtwRefB_[tableIdx]; + double TRef = watdentRefTemp_[tableIdx]; + double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]); + double cT1 = watdentCT1_[tableIdx]; + double cT2 = watdentCT2_[tableIdx]; + double Y = T[i] - TRef; + double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y); + output_b[i] = 1.0/Bw; + } + + std::fill(output_dbdp, output_dbdp + n, 0.0); + std::fill(output_dbdr, output_dbdr + n, 0.0); + } + + virtual void b(const int n, + const int* pvtRegionIdx, + const double* p, + const double* T, + const double* r, + const PhasePresence* cond, + double* output_b, + double* output_dbdp, + double* output_dbdr) const + { + if (watdentRefTemp_.empty()) { + // isothermal case + isothermalPvt_->b(n, pvtRegionIdx, p, T, r, cond, output_b, output_dbdp, output_dbdr); + return; + } + + // This changes pressure dependence of the water density, but it seems to be + // what Eclipse does. See the documentation for the WATDENT keyword in the + // Eclipse RM + for (int i = 0; i < n; ++i) { + int tableIdx = getTableIndex_(pvtRegionIdx, i); + double BwRef = pvtwRefB_[tableIdx]; + double TRef = watdentRefTemp_[tableIdx]; + double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]); + double cT1 = watdentCT1_[tableIdx]; + double cT2 = watdentCT2_[tableIdx]; + double Y = T[i] - TRef; + double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y); + output_b[i] = 1.0/Bw; + } + + std::fill(output_dbdp, output_dbdp + n, 0.0); + std::fill(output_dbdr, output_dbdr + n, 0.0); + + } + + virtual void rsSat(const int n, + const int* pvtRegionIdx, + const double* p, + double* output_rsSat, + double* output_drsSatdp) const + { + isothermalPvt_->rsSat(n, pvtRegionIdx, p, output_rsSat, output_drsSatdp); + } + + virtual void rvSat(const int n, + const int* pvtRegionIdx, + const double* p, + double* output_rvSat, + double* output_drvSatdp) const + { + isothermalPvt_->rvSat(n, pvtRegionIdx, p, output_rvSat, output_drvSatdp); + } + + virtual void R(const int n, + const int* pvtRegionIdx, + const double* p, + const double* z, + double* output_R) const + { + isothermalPvt_->R(n, pvtRegionIdx, p, z, output_R); + } + + virtual void dRdp(const int n, + const int* pvtRegionIdx, + const double* p, + const double* z, + double* output_R, + double* output_dRdp) const + { + isothermalPvt_->dRdp(n, pvtRegionIdx, p, z, output_R, output_dRdp); + } + + private: + int getTableIndex_(const int* pvtTableIdx, int cellIdx) const + { + if (!pvtTableIdx) + return 0; + return pvtTableIdx[cellIdx]; + } + + // the PVT propertied for the isothermal case + std::shared_ptr isothermalPvt_; + + // The PVT properties needed for temperature dependence. We need to store one + // value per PVT region. + std::vector viscrefPress_; + + std::vector watdentRefTemp_; + std::vector watdentCT1_; + std::vector watdentCT2_; + + std::vector pvtwRefPress_; + std::vector pvtwRefB_; + std::vector pvtwCompressibility_; + std::vector pvtwViscosity_; + std::vector pvtwViscosibility_; + + const std::vector* watvisctTables_; + }; + +} + +#endif + From 0b85af1bb6c4c0f23eceff151e409594416e6d23 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 17 Mar 2015 12:40:09 +0100 Subject: [PATCH 3/4] add a wrapper for thermal oil PVT objects --- CMakeLists_files.cmake | 1 + opm/core/props/pvt/ThermalOilPvtWrapper.hpp | 383 ++++++++++++++++++++ 2 files changed, 384 insertions(+) create mode 100644 opm/core/props/pvt/ThermalOilPvtWrapper.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 5a4d37a5..9d6e3ce5 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -357,6 +357,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/props/pvt/PvtLiveGas.hpp opm/core/props/pvt/PvtLiveOil.hpp opm/core/props/pvt/ThermalWaterPvtWrapper.hpp + opm/core/props/pvt/ThermalOilPvtWrapper.hpp opm/core/props/rock/RockBasic.hpp opm/core/props/rock/RockCompressibility.hpp opm/core/props/rock/RockFromDeck.hpp diff --git a/opm/core/props/pvt/ThermalOilPvtWrapper.hpp b/opm/core/props/pvt/ThermalOilPvtWrapper.hpp new file mode 100644 index 00000000..1ebd76c5 --- /dev/null +++ b/opm/core/props/pvt/ThermalOilPvtWrapper.hpp @@ -0,0 +1,383 @@ +/* + Copyright 2015 Andreas Lauser + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_THERMAL_OIL_PVT_WRAPPER_HPP +#define OPM_THERMAL_OIL_PVT_WRAPPER_HPP + +#include +#include + +#include + +#include + +namespace Opm +{ + /// Class which wraps another (i.e., isothermal) PVT object into one which adds + /// temperature dependence of oil + class ThermalOilPvtWrapper : public PvtInterface + { + public: + ThermalOilPvtWrapper() + {} + + + /// set the tables which specify the temperature dependence of the oil viscosity + void initFromDeck(std::shared_ptr isothermalPvt, + Opm::DeckConstPtr deck, + Opm::EclipseStateConstPtr eclipseState) + { + isothermalPvt_ = isothermalPvt; + + int numRegions; + if (deck->hasKeyword("PVTO")) + numRegions = eclipseState->getPvtoTables().size(); + else if (deck->hasKeyword("PVDO")) + numRegions = eclipseState->getPvdoTables().size(); + else if (deck->hasKeyword("PVCDO")) + numRegions = deck->getKeyword("PVCDO")->size(); + else + OPM_THROW(std::runtime_error, "Oil phase was not initialized using a known way"); + + // viscosity + if (deck->hasKeyword("VISCREF")) { + oilvisctTables_ = &eclipseState->getOilvisctTables(); + Opm::DeckKeywordConstPtr viscrefKeyword = deck->getKeyword("VISCREF"); + + assert(oilvisctTables_->size() == numRegions); + assert(viscrefKeyword->size() == numRegions); + + viscrefPress_.resize(numRegions); + viscrefRs_.resize(numRegions); + muRef_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + DeckRecordConstPtr viscrefRecord = viscrefKeyword->getRecord(regionIdx); + viscrefPress_[regionIdx] = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); + viscrefRs_[regionIdx] = viscrefRecord->getItem("REFERENCE_RS")->getSIDouble(0); + + // temperature used to calculate the reference viscosity [K]. the + // value does not really matter if the underlying PVT object really + // is isothermal... + double Tref = 273.15 + 20; + + // compute the reference viscosity using the isothermal PVT object. + double tmp1, tmp2; + isothermalPvt_->mu(1, + ®ionIdx, + &viscrefPress_[regionIdx], + &Tref, + &viscrefRs_[regionIdx], + &muRef_[regionIdx], + &tmp1, + &tmp2); + } + } + + // quantities required for density. note that we just always use the values + // for the first EOS. (since EOS != PVT region.) + tref_ = 0.0; + if (deck->hasKeyword("THERMEX1")) { + oilCompIdx_ = deck->getKeyword("OCOMPIDX")->getRecord(0)->getItem("OIL_COMPONENT_INDEX")->getInt(0) - 1; + + // always use the values of the first EOS + tref_ = deck->getKeyword("TREF")->getRecord(0)->getItem("TEMPERATURE")->getSIDouble(oilCompIdx_); + pref_ = deck->getKeyword("PREF")->getRecord(0)->getItem("PRESSURE")->getSIDouble(oilCompIdx_); + cref_ = deck->getKeyword("CREF")->getRecord(0)->getItem("COMPRESSIBILITY")->getSIDouble(oilCompIdx_); + thermex1_ = deck->getKeyword("THERMEX1")->getRecord(0)->getItem("EXPANSION_COEFF")->getSIDouble(oilCompIdx_); + } + } + + virtual void mu(const int n, + const int* pvtRegionIdx, + const double* p, + const double* T, + const double* z, + double* output_mu) const + { + if (oilvisctTables_) + // TODO: temperature dependence for viscosity depending on z + OPM_THROW(std::runtime_error, + "temperature dependent viscosity as a function of z " + "is not yet implemented!"); + + // compute the isothermal viscosity + isothermalPvt_->mu(n, pvtRegionIdx, p, T, z, output_mu); + } + + 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 + { + // compute the isothermal viscosity and its derivatives + isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, output_mu, output_dmudp, output_dmudr); + + if (!oilvisctTables_) + // isothermal case + return; + + // temperature dependence + for (int i = 0; i < n; ++i) { + int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i); + + // calculate the viscosity of the isothermal keyword for the reference + // pressure given by the VISCREF keyword. + double muRef = muRef_[regionIdx]; + + // compute the viscosity deviation due to temperature + double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]); + double alpha = muOilvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + // TODO (?): derivative of viscosity w.r.t. temperature. + } + } + + 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 + { + // compute the isothermal viscosity and its derivatives + isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, cond, output_mu, output_dmudp, output_dmudr); + + if (!oilvisctTables_) + // isothermal case + return; + + // temperature dependence + for (int i = 0; i < n; ++i) { + int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i); + + // calculate the viscosity of the isothermal keyword for the reference + // pressure given by the VISCREF keyword. + double muRef = muRef_[regionIdx]; + + // compute the viscosity deviation due to temperature + double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]); + double alpha = muOilvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + // TODO (?): derivative of viscosity w.r.t. temperature. + } + } + + virtual void B(const int n, + const int* pvtRegionIdx, + const double* p, + const double* T, + const double* z, + double* output_B) const + { + // isothermal case + isothermalPvt_->B(n, pvtRegionIdx, p, T, z, output_B); + + if (thermex1_ <= 0.0) + // isothermal case + return; + + // deal with the temperature dependence of the oil phase. we use equation + // (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref + // using the isothermal keyword instead of using the value for the + // components, so the oil compressibility is already dealt with there. Note + // that we only do the part for the oil component here, the part for + // dissolved gas is ignored so far. + double cT1 = thermex1_; + double TRef = tref_; + for (int i = 0; i < n; ++i) { + double alpha = (1 + cT1*(T[i] - TRef)); + output_B[i] *= alpha; + } + } + + 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 + { + isothermalPvt_->dBdp(n, pvtRegionIdx, p, T, z, output_B, output_dBdp); + + if (thermex1_ <= 0.0) + // isothermal case + return; + + // deal with the temperature dependence of the oil phase. we use equation + // (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref + // using the isothermal keyword instead of using the value for the + // components, so the oil compressibility is already dealt with there. Note + // that we only do the part for the oil component here, the part for + // dissolved gas is ignored so far. + double cT1 = thermex1_; + double TRef = tref_; + for (int i = 0; i < n; ++i) { + double alpha = (1 + cT1*(T[i] - TRef)); + output_B[i] *= alpha; + output_dBdp[i] *= alpha; + } + } + + 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 + { + isothermalPvt_->b(n, pvtRegionIdx, p, T, r, output_b, output_dbdp, output_dbdr); + + if (thermex1_ <= 0.0) + // isothermal case + return; + + // deal with the temperature dependence of the oil phase. we use equation + // (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref + // using the isothermal keyword instead of using the value for the + // components, so the oil compressibility is already dealt with there. Note + // that we only do the part for the oil component here, the part for + // dissolved gas is ignored so far. + double cT1 = thermex1_; + double TRef = tref_; + for (int i = 0; i < n; ++i) { + double alpha = 1.0/(1 + cT1*(T[i] - TRef)); + output_b[i] *= alpha; + output_dbdp[i] *= alpha; + output_dbdr[i] *= alpha; + } + } + + virtual void b(const int n, + const int* pvtRegionIdx, + const double* p, + const double* T, + const double* r, + const PhasePresence* cond, + double* output_b, + double* output_dbdp, + double* output_dbdr) const + { + isothermalPvt_->b(n, pvtRegionIdx, p, T, r, cond, output_b, output_dbdp, output_dbdr); + + if (thermex1_ <= 0.0) + // isothermal case + return; + + // deal with the temperature dependence of the oil phase. we use equation + // (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref + // using the isothermal keyword instead of using the value for the + // components, so the oil compressibility is already dealt with there. Note + // that we only do the part for the oil component here, the part for + // dissolved gas is ignored so far. + double cT1 = thermex1_; + double TRef = tref_; + for (int i = 0; i < n; ++i) { + double alpha = 1.0/(1 + cT1*(T[i] - TRef)); + output_b[i] *= alpha; + output_dbdp[i] *= alpha; + output_dbdr[i] *= alpha; + } + } + + virtual void rsSat(const int n, + const int* pvtRegionIdx, + const double* p, + double* output_rsSat, + double* output_drsSatdp) const + { + isothermalPvt_->rsSat(n, pvtRegionIdx, p, output_rsSat, output_drsSatdp); + } + + virtual void rvSat(const int n, + const int* pvtRegionIdx, + const double* p, + double* output_rvSat, + double* output_drvSatdp) const + { + isothermalPvt_->rvSat(n, pvtRegionIdx, p, output_rvSat, output_drvSatdp); + } + + virtual void R(const int n, + const int* pvtRegionIdx, + const double* p, + const double* z, + double* output_R) const + { + isothermalPvt_->R(n, pvtRegionIdx, p, z, output_R); + } + + virtual void dRdp(const int n, + const int* pvtRegionIdx, + const double* p, + const double* z, + double* output_R, + double* output_dRdp) const + { + isothermalPvt_->dRdp(n, pvtRegionIdx, p, z, output_R, output_dRdp); + } + + private: + int getPvtRegionIndex_(const int* pvtRegionIdx, int cellIdx) const + { + if (!pvtRegionIdx) + return 0; + return pvtRegionIdx[cellIdx]; + } + + // the PVT propertied for the isothermal case + std::shared_ptr isothermalPvt_; + + // The PVT properties needed for temperature dependence of the viscosity. We need + // to store one value per PVT region. + std::vector viscrefPress_; + std::vector viscrefRs_; + std::vector muRef_; + + const std::vector* oilvisctTables_; + + // The PVT properties needed for temperature dependence of the density. This is + // specified as one value per EOS in the manual, but we unconditionally use the + // expansion coefficient of the first EOS... + int oilCompIdx_; + double tref_; + double pref_; + double cref_; + double thermex1_; + }; + +} + +#endif + From a386621495fd9378dedde917b3517f27f9342991 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 17 Mar 2015 12:40:14 +0100 Subject: [PATCH 4/4] add a wrapper for thermal gas PVT objects --- CMakeLists_files.cmake | 1 + opm/core/props/pvt/ThermalGasPvtWrapper.hpp | 310 ++++++++++++++++++++ 2 files changed, 311 insertions(+) create mode 100644 opm/core/props/pvt/ThermalGasPvtWrapper.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 9d6e3ce5..14cc1aef 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -358,6 +358,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/props/pvt/PvtLiveOil.hpp opm/core/props/pvt/ThermalWaterPvtWrapper.hpp opm/core/props/pvt/ThermalOilPvtWrapper.hpp + opm/core/props/pvt/ThermalGasPvtWrapper.hpp opm/core/props/rock/RockBasic.hpp opm/core/props/rock/RockCompressibility.hpp opm/core/props/rock/RockFromDeck.hpp diff --git a/opm/core/props/pvt/ThermalGasPvtWrapper.hpp b/opm/core/props/pvt/ThermalGasPvtWrapper.hpp new file mode 100644 index 00000000..39b541bb --- /dev/null +++ b/opm/core/props/pvt/ThermalGasPvtWrapper.hpp @@ -0,0 +1,310 @@ +/* + Copyright 2015 Andreas Lauser + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef OPM_THERMAL_GAS_PVT_WRAPPER_HPP +#define OPM_THERMAL_GAS_PVT_WRAPPER_HPP + +#include +#include + +#include + +#include + +namespace Opm +{ + /// Class which wraps another (i.e., isothermal) PVT object into one which adds + /// temperature dependence of gas + class ThermalGasPvtWrapper : public PvtInterface + { + public: + ThermalGasPvtWrapper() + {} + + + /// extract the quantities needed specify the temperature dependence of the gas + /// viscosity and density from the deck + void initFromDeck(std::shared_ptr isothermalPvt, + Opm::DeckConstPtr deck, + Opm::EclipseStateConstPtr eclipseState) + { + isothermalPvt_ = isothermalPvt; + + int numRegions; + if (deck->hasKeyword("PVTG")) + numRegions = eclipseState->getPvtgTables().size(); + else if (deck->hasKeyword("PVDG")) + numRegions = eclipseState->getPvdgTables().size(); + else + OPM_THROW(std::runtime_error, "Gas phase was not initialized using a known way"); + + // viscosity + if (deck->hasKeyword("GASVISCT")) { + gasvisctTables_ = &eclipseState->getGasvisctTables(); + assert(gasvisctTables_->size() == numRegions); + + gasCompIdx_ = deck->getKeyword("GCOMPIDX")->getRecord(0)->getItem("GAS_COMPONENT_INDEX")->getInt(0) - 1; + gasvisctColumnName_ = "Viscosity"+std::to_string(static_cast(gasCompIdx_)); + } + + // density + if (deck->hasKeyword("TREF")) { + tref_ = deck->getKeyword("TREF")->getRecord(0)->getItem("TEMPERATURE")->getSIDouble(0); + } + } + + virtual void mu(const int n, + const int* pvtRegionIdx, + const double* p, + const double* T, + const double* z, + double* output_mu) const + { + if (gasvisctTables_) + // TODO: temperature dependence for viscosity depending on z + OPM_THROW(std::runtime_error, + "temperature dependent viscosity as a function of z " + "is not yet implemented!"); + + // compute the isothermal viscosity + isothermalPvt_->mu(n, pvtRegionIdx, p, T, z, output_mu); + } + + 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 + { + if (gasvisctTables_ != 0) { + for (int i = 0; i < n; ++i) { + // temperature dependence of the gas phase. this assumes that the gas + // component index has been set properly, and it also looses the + // pressure dependence of gas. (This does not make much sense, but it + // seems to be what the documentation for the GASVISCT keyword in the + // RM says.) + int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i); + double muGasvisct = (*gasvisctTables_)[regionIdx].evaluate(gasvisctColumnName_, T[i]); + + output_mu[i] = muGasvisct; + output_dmudp[i] = 0.0; + output_dmudr[i] = 0.0; + + // TODO (?): derivative of gas viscosity w.r.t. temperature. + } + } + else { + // compute the isothermal viscosity and its derivatives + isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, output_mu, output_dmudp, output_dmudr); + } + } + + 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 + { + if (gasvisctTables_ != 0) { + for (int i = 0; i < n; ++i) { + // temperature dependence of the gas phase. this assumes that the gas + // component index has been set properly, and it also looses the + // pressure dependence of gas. (This does not make much sense, but it + // seems to be what the documentation for the GASVISCT keyword in the + // RM says.) + int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i); + double muGasvisct = (*gasvisctTables_)[regionIdx].evaluate(gasvisctColumnName_, T[i]); + + output_mu[i] = muGasvisct; + output_dmudp[i] = 0.0; + output_dmudr[i] = 0.0; + + // TODO (?): derivative of gas viscosity w.r.t. temperature. + } + } + else { + // compute the isothermal viscosity and its derivatives + isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, cond, output_mu, output_dmudp, output_dmudr); + } + } + + virtual void B(const int n, + const int* pvtRegionIdx, + const double* p, + const double* T, + const double* z, + double* output_B) const + { + // isothermal case + isothermalPvt_->B(n, pvtRegionIdx, p, T, z, output_B); + + if (tref_ > 0.0) { + // the Eclipse TD/RM do not explicitly specify the relation of the gas + // density and the temperature, but equation (69.49) (for Eclipse 2011.1) + // implies that the temperature dependence of the gas phase is rho(T, p) = + // rho(tref_, p)/tref_*T ... + for (int i = 0; i < n; ++i) { + double alpha = tref_/T[i]; + output_B[i] *= alpha; + } + } + } + + 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 + { + isothermalPvt_->dBdp(n, pvtRegionIdx, p, T, z, output_B, output_dBdp); + + if (tref_ > 0.0) { + // the Eclipse TD/RM do not explicitly specify the relation of the gas + // density and the temperature, but equation (69.49) (for Eclipse 2011.1) + // implies that the temperature dependence of the gas phase is rho(T, p) = + // rho(tref_, p)/tref_*T ... + for (int i = 0; i < n; ++i) { + double alpha = tref_/T[i]; + output_B[i] *= alpha; + output_dBdp[i] *= alpha; + } + } + } + + 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 + { + isothermalPvt_->b(n, pvtRegionIdx, p, T, r, output_b, output_dbdp, output_dbdr); + + if (tref_ > 0.0) { + // the Eclipse TD/RM do not explicitly specify the relation of the gas + // density and the temperature, but equation (69.49) (for Eclipse 2011.1) + // implies that the temperature dependence of the gas phase is rho(T, p) = + // rho(tref_, p)/tref_*T ... + for (int i = 0; i < n; ++i) { + double alpha = T[i]/tref_; + output_b[i] *= alpha; + output_dbdp[i] *= alpha; + output_dbdr[i] *= alpha; + } + } + } + + virtual void b(const int n, + const int* pvtRegionIdx, + const double* p, + const double* T, + const double* r, + const PhasePresence* cond, + double* output_b, + double* output_dbdp, + double* output_dbdr) const + { + isothermalPvt_->b(n, pvtRegionIdx, p, T, r, cond, output_b, output_dbdp, output_dbdr); + + if (tref_ > 0.0) { + // the Eclipse TD/RM do not explicitly specify the relation of the gas + // density and the temperature, but equation (69.49) (for Eclipse 2011.1) + // implies that the temperature dependence of the gas phase is rho(T, p) = + // rho(tref_, p)/tref_*T ... + for (int i = 0; i < n; ++i) { + double alpha = T[i]/tref_; + output_b[i] *= alpha; + output_dbdp[i] *= alpha; + output_dbdr[i] *= alpha; + } + } + } + + virtual void rsSat(const int n, + const int* pvtRegionIdx, + const double* p, + double* output_rsSat, + double* output_drsSatdp) const + { + isothermalPvt_->rsSat(n, pvtRegionIdx, p, output_rsSat, output_drsSatdp); + } + + virtual void rvSat(const int n, + const int* pvtRegionIdx, + const double* p, + double* output_rvSat, + double* output_drvSatdp) const + { + isothermalPvt_->rvSat(n, pvtRegionIdx, p, output_rvSat, output_drvSatdp); + } + + virtual void R(const int n, + const int* pvtRegionIdx, + const double* p, + const double* z, + double* output_R) const + { + isothermalPvt_->R(n, pvtRegionIdx, p, z, output_R); + } + + virtual void dRdp(const int n, + const int* pvtRegionIdx, + const double* p, + const double* z, + double* output_R, + double* output_dRdp) const + { + isothermalPvt_->dRdp(n, pvtRegionIdx, p, z, output_R, output_dRdp); + } + + private: + int getPvtRegionIndex_(const int* pvtRegionIdx, int cellIdx) const + { + if (!pvtRegionIdx) + return 0; + return pvtRegionIdx[cellIdx]; + } + + // the PVT propertied for the isothermal case + std::shared_ptr isothermalPvt_; + + // The PVT properties needed for temperature dependence of the viscosity. We need + // to store one value per PVT region. + const std::vector* gasvisctTables_; + std::string gasvisctColumnName_; + int gasCompIdx_; + + // The PVT properties needed for temperature dependence of the density. + double tref_; + }; + +} + +#endif +