20#ifndef OPM_THERMAL_WATER_PVT_WRAPPER_HPP
21#define OPM_THERMAL_WATER_PVT_WRAPPER_HPP
23#include <opm/core/props/pvt/PvtInterface.hpp>
24#include <opm/common/ErrorMacros.hpp>
26#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
27#include <opm/parser/eclipse/EclipseState/Tables/WatvisctTable.hpp>
43 void initFromDeck(std::shared_ptr<const PvtInterface> isothermalPvt,
44 const Opm::Deck& deck,
45 const Opm::EclipseState& eclipseState)
47 isothermalPvt_ = isothermalPvt;
51 const auto& pvtwKeyword = deck->getKeyword(
"PVTW");
52 int numRegions = pvtwKeyword.size();
53 pvtwRefPress_.resize(numRegions);
54 pvtwRefB_.resize(numRegions);
55 pvtwCompressibility_.resize(numRegions);
56 pvtwViscosity_.resize(numRegions);
57 pvtwViscosibility_.resize(numRegions);
58 for (
int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
59 const auto& pvtwRecord = pvtwKeyword.getRecord(regionIdx);
60 pvtwRefPress_[regionIdx] = pvtwRecord.getItem(
"P_REF").getSIDouble(0);
61 pvtwRefB_[regionIdx] = pvtwRecord.getItem(
"WATER_VOL_FACTOR").getSIDouble(0);
62 pvtwViscosity_[regionIdx] = pvtwRecord.getItem(
"WATER_VISCOSITY").getSIDouble(0);
63 pvtwViscosibility_[regionIdx] = pvtwRecord.getItem(
"WATER_VISCOSIBILITY").getSIDouble(0);
68 if (deck->hasKeyword(
"VISCREF")) {
69 auto tables = eclipseState->getTableManager();
70 watvisctTables_ = &tables->getWatvisctTables();
71 const auto& viscrefKeyword = deck->getKeyword(
"VISCREF");
73 assert(
int(watvisctTables_->size()) == numRegions);
74 assert(
int(viscrefKeyword.size()) == numRegions);
76 viscrefPress_.resize(numRegions);
77 for (
int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
78 const auto& viscrefRecord = viscrefKeyword.getRecord(regionIdx);
80 viscrefPress_[regionIdx] = viscrefRecord.getItem(
"REFERENCE_PRESSURE").getSIDouble(0);
85 if (deck->hasKeyword(
"WATDENT")) {
86 const auto& watdentKeyword = deck->getKeyword(
"WATDENT");
88 assert(
int(watdentKeyword.size()) == numRegions);
90 watdentRefTemp_.resize(numRegions);
91 watdentCT1_.resize(numRegions);
92 watdentCT2_.resize(numRegions);
93 for (
int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
94 const auto& watdentRecord = watdentKeyword.getRecord(regionIdx);
96 watdentRefTemp_[regionIdx] = watdentRecord.getItem(
"REFERENCE_TEMPERATURE").getSIDouble(0);
97 watdentCT1_[regionIdx] = watdentRecord.getItem(
"EXPANSION_COEFF_LINEAR").getSIDouble(0);
98 watdentCT2_[regionIdx] = watdentRecord.getItem(
"EXPANSION_COEFF_QUADRATIC").getSIDouble(0);
103 virtual void mu(
const int n,
104 const int* pvtRegionIdx,
108 double* output_mu)
const
112 OPM_THROW(std::runtime_error,
113 "temperature dependent viscosity as a function of z "
114 "is not yet implemented!");
117 isothermalPvt_->mu(n, pvtRegionIdx, p, T, z, output_mu);
120 virtual void mu(
const int n,
121 const int* pvtRegionIdx,
126 double* output_dmudp,
127 double* output_dmudr)
const
130 isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, output_mu, output_dmudp, output_dmudr);
132 if (!watvisctTables_)
137 for (
int i = 0; i < n; ++i) {
138 int tableIdx = getTableIndex_(pvtRegionIdx, i);
142 double x = -pvtwViscosibility_[tableIdx]*(viscrefPress_[tableIdx] - pvtwRefPress_[tableIdx]);
143 double muRef = pvtwViscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
148 const WatvisctTable& watVisctTable = watvisctTables_->getTable<WatvisctTable>(tableIdx);
149 double muWatvisct = watVisctTable.evaluate(
"Viscosity", T[i]);
150 alpha = muWatvisct/muRef;
153 output_mu[i] *= alpha;
154 output_dmudp[i] *= alpha;
155 output_dmudr[i] *= alpha;
160 virtual void mu(
const int n,
161 const int* pvtRegionIdx,
167 double* output_dmudp,
168 double* output_dmudr)
const
171 isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, cond, output_mu, output_dmudp, output_dmudr);
173 if (!watvisctTables_)
178 for (
int i = 0; i < n; ++i) {
179 int tableIdx = getTableIndex_(pvtRegionIdx, i);
183 double x = -pvtwViscosibility_[tableIdx]*(viscrefPress_[tableIdx] - pvtwRefPress_[tableIdx]);
184 double muRef = pvtwViscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
189 const WatvisctTable& watVisctTable = watvisctTables_->getTable<WatvisctTable>(tableIdx);
190 double muWatvisct = watVisctTable.evaluate(
"Viscosity", T[i]);
191 alpha = muWatvisct/muRef;
193 output_mu[i] *= alpha;
194 output_dmudp[i] *= alpha;
195 output_dmudr[i] *= alpha;
200 virtual void B(
const int n,
201 const int* pvtRegionIdx,
205 double* output_B)
const
207 if (watdentRefTemp_.empty()) {
209 isothermalPvt_->B(n, pvtRegionIdx, p, T, z, output_B);
216 for (
int i = 0; i < n; ++i) {
217 int tableIdx = getTableIndex_(pvtRegionIdx, i);
218 double BwRef = pvtwRefB_[tableIdx];
219 double TRef = watdentRefTemp_[tableIdx];
220 double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
221 double cT1 = watdentCT1_[tableIdx];
222 double cT2 = watdentCT2_[tableIdx];
223 double Y = T[i] - TRef;
224 double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
230 const int* pvtRegionIdx,
235 double* output_dBdp)
const
237 if (watdentRefTemp_.empty()) {
239 isothermalPvt_->dBdp(n, pvtRegionIdx, p, T, z, output_B, output_dBdp);
246 for (
int i = 0; i < n; ++i) {
247 int tableIdx = getTableIndex_(pvtRegionIdx, i);
248 double BwRef = pvtwRefB_[tableIdx];
249 double TRef = watdentRefTemp_[tableIdx];
250 double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
251 double cT1 = watdentCT1_[tableIdx];
252 double cT2 = watdentCT2_[tableIdx];
253 double Y = T[i] - TRef;
254 double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
258 std::fill(output_dBdp, output_dBdp + n, 0.0);
261 virtual void b(
const int n,
262 const int* pvtRegionIdx,
268 double* output_dbdr)
const
270 if (watdentRefTemp_.empty()) {
272 isothermalPvt_->b(n, pvtRegionIdx, p, T, r, output_b, output_dbdp, output_dbdr);
279 for (
int i = 0; i < n; ++i) {
280 int tableIdx = getTableIndex_(pvtRegionIdx, i);
281 double BwRef = pvtwRefB_[tableIdx];
282 double TRef = watdentRefTemp_[tableIdx];
283 double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
284 double cT1 = watdentCT1_[tableIdx];
285 double cT2 = watdentCT2_[tableIdx];
286 double Y = T[i] - TRef;
287 double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
288 output_b[i] = 1.0/Bw;
291 std::fill(output_dbdp, output_dbdp + n, 0.0);
292 std::fill(output_dbdr, output_dbdr + n, 0.0);
295 virtual void b(
const int n,
296 const int* pvtRegionIdx,
303 double* output_dbdr)
const
305 if (watdentRefTemp_.empty()) {
307 isothermalPvt_->b(n, pvtRegionIdx, p, T, r, cond, output_b, output_dbdp, output_dbdr);
314 for (
int i = 0; i < n; ++i) {
315 int tableIdx = getTableIndex_(pvtRegionIdx, i);
316 double BwRef = pvtwRefB_[tableIdx];
317 double TRef = watdentRefTemp_[tableIdx];
318 double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
319 double cT1 = watdentCT1_[tableIdx];
320 double cT2 = watdentCT2_[tableIdx];
321 double Y = T[i] - TRef;
322 double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
323 output_b[i] = 1.0/Bw;
326 std::fill(output_dbdp, output_dbdp + n, 0.0);
327 std::fill(output_dbdr, output_dbdr + n, 0.0);
332 const int* pvtRegionIdx,
334 double* output_rsSat,
335 double* output_drsSatdp)
const
337 isothermalPvt_->rsSat(n, pvtRegionIdx, p, output_rsSat, output_drsSatdp);
341 const int* pvtRegionIdx,
343 double* output_rvSat,
344 double* output_drvSatdp)
const
346 isothermalPvt_->rvSat(n, pvtRegionIdx, p, output_rvSat, output_drvSatdp);
349 virtual void R(
const int n,
350 const int* pvtRegionIdx,
353 double* output_R)
const
355 isothermalPvt_->R(n, pvtRegionIdx, p, z, output_R);
359 const int* pvtRegionIdx,
363 double* output_dRdp)
const
365 isothermalPvt_->dRdp(n, pvtRegionIdx, p, z, output_R, output_dRdp);
369 int getTableIndex_(
const int* pvtTableIdx,
int cellIdx)
const
373 return pvtTableIdx[cellIdx];
377 std::shared_ptr<const PvtInterface> isothermalPvt_;
381 std::vector<double> viscrefPress_;
383 std::vector<double> watdentRefTemp_;
384 std::vector<double> watdentCT1_;
385 std::vector<double> watdentCT2_;
387 std::vector<double> pvtwRefPress_;
388 std::vector<double> pvtwRefB_;
389 std::vector<double> pvtwCompressibility_;
390 std::vector<double> pvtwViscosity_;
391 std::vector<double> pvtwViscosibility_;
393 const TableContainer* watvisctTables_;
Definition: BlackoilPhases.hpp:59
Definition: ThermalWaterPvtWrapper.hpp:36
void initFromDeck(std::shared_ptr< const PvtInterface > isothermalPvt, const Opm::Deck &deck, const Opm::EclipseState &eclipseState)
set the tables which specify the temperature dependence of the water viscosity
Definition: ThermalWaterPvtWrapper.hpp:43
virtual void dRdp(const int n, const int *pvtRegionIdx, const double *p, const double *z, double *output_R, double *output_dRdp) const
Definition: ThermalWaterPvtWrapper.hpp:358
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
Definition: ThermalWaterPvtWrapper.hpp:295
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
Definition: ThermalWaterPvtWrapper.hpp:229
virtual void R(const int n, const int *pvtRegionIdx, const double *p, const double *z, double *output_R) const
Definition: ThermalWaterPvtWrapper.hpp:349
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
Definition: ThermalWaterPvtWrapper.hpp:261
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
Definition: ThermalWaterPvtWrapper.hpp:120
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
Definition: ThermalWaterPvtWrapper.hpp:160
virtual void rsSat(const int n, const int *pvtRegionIdx, const double *p, double *output_rsSat, double *output_drsSatdp) const
Definition: ThermalWaterPvtWrapper.hpp:331
ThermalWaterPvtWrapper()
Definition: ThermalWaterPvtWrapper.hpp:38
virtual void mu(const int n, const int *pvtRegionIdx, const double *p, const double *T, const double *z, double *output_mu) const
Definition: ThermalWaterPvtWrapper.hpp:103
virtual void rvSat(const int n, const int *pvtRegionIdx, const double *p, double *output_rvSat, double *output_drvSatdp) const
Definition: ThermalWaterPvtWrapper.hpp:340
virtual void B(const int n, const int *pvtRegionIdx, const double *p, const double *T, const double *z, double *output_B) const
Definition: ThermalWaterPvtWrapper.hpp:200
Definition: AnisotropicEikonal.hpp:44