All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Pages
ThermalGasPvtWrapper.hpp
1 /*
2  Copyright 2015 Andreas Lauser
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 #ifndef OPM_THERMAL_GAS_PVT_WRAPPER_HPP
20 #define OPM_THERMAL_GAS_PVT_WRAPPER_HPP
21 
22 #include <opm/core/props/pvt/PvtInterface.hpp>
23 #include <opm/common/ErrorMacros.hpp>
24 
25 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
26 #include <opm/parser/eclipse/EclipseState/Tables/GasvisctTable.hpp>
27 
28 #include <vector>
29 
30 namespace Opm
31 {
34  class ThermalGasPvtWrapper : public PvtInterface
35  {
36  public:
38  {}
39 
40 
43  void initFromDeck(std::shared_ptr<const PvtInterface> isothermalPvt,
44  const Opm::Deck& deck,
45  const Opm::EclipseState& eclipseState)
46  {
47  isothermalPvt_ = isothermalPvt;
48  auto tables = eclipseState->getTableManager();
49  int numRegions;
50  if (deck->hasKeyword("PVTG"))
51  numRegions = tables->getPvtgTables().size();
52  else if (deck->hasKeyword("PVDG"))
53  numRegions = tables->getPvdgTables().size();
54  else
55  OPM_THROW(std::runtime_error, "Gas phase was not initialized using a known way");
56 
57  // viscosity
58  if (deck->hasKeyword("GASVISCT")) {
59  gasvisctTables_ = &tables->getGasvisctTables();
60  assert(int(gasvisctTables_->size()) == numRegions);
61  static_cast<void>(numRegions); //Silence compiler warning
62 
63  gasCompIdx_ = deck->getKeyword("GCOMPIDX").getRecord(0).getItem("GAS_COMPONENT_INDEX").get< int >(0) - 1;
64  gasvisctColumnName_ = "Viscosity"+std::to_string(static_cast<long long>(gasCompIdx_));
65  }
66 
67  // density
68  if (deck->hasKeyword("TREF")) {
69  tref_ = deck->getKeyword("TREF").getRecord(0).getItem("TEMPERATURE").getSIDouble(0);
70  }
71  }
72 
73  virtual void mu(const int n,
74  const int* pvtRegionIdx,
75  const double* p,
76  const double* T,
77  const double* z,
78  double* output_mu) const
79  {
80  if (gasvisctTables_)
81  // TODO: temperature dependence for viscosity depending on z
82  OPM_THROW(std::runtime_error,
83  "temperature dependent viscosity as a function of z "
84  "is not yet implemented!");
85 
86  // compute the isothermal viscosity
87  isothermalPvt_->mu(n, pvtRegionIdx, p, T, z, output_mu);
88  }
89 
90  virtual void mu(const int n,
91  const int* pvtRegionIdx,
92  const double* p,
93  const double* T,
94  const double* r,
95  double* output_mu,
96  double* output_dmudp,
97  double* output_dmudr) const
98  {
99  if (gasvisctTables_ != 0) {
100  for (int i = 0; i < n; ++i) {
101  // temperature dependence of the gas phase. this assumes that the gas
102  // component index has been set properly, and it also looses the
103  // pressure dependence of gas. (This does not make much sense, but it
104  // seems to be what the documentation for the GASVISCT keyword in the
105  // RM says.)
106 
107 
108  int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i);
109  double muGasvisct;
110  {
111  const GasvisctTable& gasvisctTable = gasvisctTables_->getTable<GasvisctTable>(regionIdx);
112  muGasvisct = gasvisctTable.evaluate(gasvisctColumnName_, T[i]);
113  }
114 
115  output_mu[i] = muGasvisct;
116  output_dmudp[i] = 0.0;
117  output_dmudr[i] = 0.0;
118 
119  // TODO (?): derivative of gas viscosity w.r.t. temperature.
120  }
121  }
122  else {
123  // compute the isothermal viscosity and its derivatives
124  isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, output_mu, output_dmudp, output_dmudr);
125  }
126  }
127 
128  virtual void mu(const int n,
129  const int* pvtRegionIdx,
130  const double* p,
131  const double* T,
132  const double* r,
133  const PhasePresence* cond,
134  double* output_mu,
135  double* output_dmudp,
136  double* output_dmudr) const
137  {
138  if (gasvisctTables_ != 0) {
139  for (int i = 0; i < n; ++i) {
140  // temperature dependence of the gas phase. this assumes that the gas
141  // component index has been set properly, and it also looses the
142  // pressure dependence of gas. (This does not make much sense, but it
143  // seems to be what the documentation for the GASVISCT keyword in the
144  // RM says.)
145  int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i);
146  double muGasvisct;
147  {
148  const GasvisctTable& gasvisctTable = gasvisctTables_->getTable<GasvisctTable>(regionIdx);
149  muGasvisct = gasvisctTable.evaluate(gasvisctColumnName_, T[i]);
150  }
151 
152  output_mu[i] = muGasvisct;
153  output_dmudp[i] = 0.0;
154  output_dmudr[i] = 0.0;
155 
156  // TODO (?): derivative of gas viscosity w.r.t. temperature.
157  }
158  }
159  else {
160  // compute the isothermal viscosity and its derivatives
161  isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, cond, output_mu, output_dmudp, output_dmudr);
162  }
163  }
164 
165  virtual void B(const int n,
166  const int* pvtRegionIdx,
167  const double* p,
168  const double* T,
169  const double* z,
170  double* output_B) const
171  {
172  // isothermal case
173  isothermalPvt_->B(n, pvtRegionIdx, p, T, z, output_B);
174 
175  if (tref_ > 0.0) {
176  // the Eclipse TD/RM do not explicitly specify the relation of the gas
177  // density and the temperature, but equation (69.49) (for Eclipse 2011.1)
178  // implies that the temperature dependence of the gas phase is rho(T, p) =
179  // rho(tref_, p)/tref_*T ...
180  for (int i = 0; i < n; ++i) {
181  double alpha = tref_/T[i];
182  output_B[i] *= alpha;
183  }
184  }
185  }
186 
187  virtual void dBdp(const int n,
188  const int* pvtRegionIdx,
189  const double* p,
190  const double* T,
191  const double* z,
192  double* output_B,
193  double* output_dBdp) const
194  {
195  isothermalPvt_->dBdp(n, pvtRegionIdx, p, T, z, output_B, output_dBdp);
196 
197  if (tref_ > 0.0) {
198  // the Eclipse TD/RM do not explicitly specify the relation of the gas
199  // density and the temperature, but equation (69.49) (for Eclipse 2011.1)
200  // implies that the temperature dependence of the gas phase is rho(T, p) =
201  // rho(tref_, p)/tref_*T ...
202  for (int i = 0; i < n; ++i) {
203  double alpha = tref_/T[i];
204  output_B[i] *= alpha;
205  output_dBdp[i] *= alpha;
206  }
207  }
208  }
209 
210  virtual void b(const int n,
211  const int* pvtRegionIdx,
212  const double* p,
213  const double* T,
214  const double* r,
215  double* output_b,
216  double* output_dbdp,
217  double* output_dbdr) const
218  {
219  isothermalPvt_->b(n, pvtRegionIdx, p, T, r, output_b, output_dbdp, output_dbdr);
220 
221  if (tref_ > 0.0) {
222  // the Eclipse TD/RM do not explicitly specify the relation of the gas
223  // density and the temperature, but equation (69.49) (for Eclipse 2011.1)
224  // implies that the temperature dependence of the gas phase is rho(T, p) =
225  // rho(tref_, p)/tref_*T ...
226  for (int i = 0; i < n; ++i) {
227  double alpha = T[i]/tref_;
228  output_b[i] *= alpha;
229  output_dbdp[i] *= alpha;
230  output_dbdr[i] *= alpha;
231  }
232  }
233  }
234 
235  virtual void b(const int n,
236  const int* pvtRegionIdx,
237  const double* p,
238  const double* T,
239  const double* r,
240  const PhasePresence* cond,
241  double* output_b,
242  double* output_dbdp,
243  double* output_dbdr) const
244  {
245  isothermalPvt_->b(n, pvtRegionIdx, p, T, r, cond, output_b, output_dbdp, output_dbdr);
246 
247  if (tref_ > 0.0) {
248  // the Eclipse TD/RM do not explicitly specify the relation of the gas
249  // density and the temperature, but equation (69.49) (for Eclipse 2011.1)
250  // implies that the temperature dependence of the gas phase is rho(T, p) =
251  // rho(tref_, p)/tref_*T ...
252  for (int i = 0; i < n; ++i) {
253  double alpha = T[i]/tref_;
254  output_b[i] *= alpha;
255  output_dbdp[i] *= alpha;
256  output_dbdr[i] *= alpha;
257  }
258  }
259  }
260 
261  virtual void rsSat(const int n,
262  const int* pvtRegionIdx,
263  const double* p,
264  double* output_rsSat,
265  double* output_drsSatdp) const
266  {
267  isothermalPvt_->rsSat(n, pvtRegionIdx, p, output_rsSat, output_drsSatdp);
268  }
269 
270  virtual void rvSat(const int n,
271  const int* pvtRegionIdx,
272  const double* p,
273  double* output_rvSat,
274  double* output_drvSatdp) const
275  {
276  isothermalPvt_->rvSat(n, pvtRegionIdx, p, output_rvSat, output_drvSatdp);
277  }
278 
279  virtual void R(const int n,
280  const int* pvtRegionIdx,
281  const double* p,
282  const double* z,
283  double* output_R) const
284  {
285  isothermalPvt_->R(n, pvtRegionIdx, p, z, output_R);
286  }
287 
288  virtual void dRdp(const int n,
289  const int* pvtRegionIdx,
290  const double* p,
291  const double* z,
292  double* output_R,
293  double* output_dRdp) const
294  {
295  isothermalPvt_->dRdp(n, pvtRegionIdx, p, z, output_R, output_dRdp);
296  }
297 
298  private:
299  int getPvtRegionIndex_(const int* pvtRegionIdx, int cellIdx) const
300  {
301  if (!pvtRegionIdx)
302  return 0;
303  return pvtRegionIdx[cellIdx];
304  }
305 
306  // the PVT propertied for the isothermal case
307  std::shared_ptr<const PvtInterface> isothermalPvt_;
308 
309  // The PVT properties needed for temperature dependence of the viscosity. We need
310  // to store one value per PVT region.
311  const TableContainer* gasvisctTables_;
312  std::string gasvisctColumnName_;
313  int gasCompIdx_;
314 
315  // The PVT properties needed for temperature dependence of the density.
316  double tref_;
317  };
318 
319 }
320 
321 #endif
322 
void initFromDeck(std::shared_ptr< const PvtInterface > isothermalPvt, const Opm::Deck &deck, const Opm::EclipseState &eclipseState)
extract the quantities needed specify the temperature dependence of the gas viscosity and density fro...
Definition: ThermalGasPvtWrapper.hpp:43
Class which wraps another (i.e., isothermal) PVT object into one which adds temperature dependence of...
Definition: ThermalGasPvtWrapper.hpp:34