My Project
EclEpsScalingPoints.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
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 2 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  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_ECL_EPS_SCALING_POINTS_HPP
28 #define OPM_ECL_EPS_SCALING_POINTS_HPP
29 
30 #include "EclEpsConfig.hpp"
31 #include "EclEpsGridProperties.hpp"
32 
33 #if HAVE_ECL_INPUT
34 #include <opm/input/eclipse/Deck/Deck.hpp>
35 #include <opm/input/eclipse/Deck/DeckRecord.hpp>
36 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
37 #include <opm/input/eclipse/EclipseState/Runspec.hpp>
38 #include <opm/input/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.hpp>
39 #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
40 #endif
41 
43 
44 #include <array>
45 #include <vector>
46 #include <string>
47 #include <iostream>
48 #include <cassert>
49 #include <algorithm>
50 
51 namespace Opm {
52 
60 template <class Scalar>
62 {
63  // connate saturations
64  Scalar Swl; // water
65  Scalar Sgl; // gas
66 
67  // critical saturations
68  Scalar Swcr; // water
69  Scalar Sgcr; // gas
70  Scalar Sowcr; // oil for the oil-water system
71  Scalar Sogcr; // oil for the gas-oil system
72 
73  // maximum saturations
74  Scalar Swu; // water
75  Scalar Sgu; // gas
76 
77  // maximum capillary pressures
78  Scalar maxPcow; // maximum capillary pressure of the oil-water system
79  Scalar maxPcgo; // maximum capillary pressure of the gas-oil system
80 
81  // the Leverett capillary pressure scaling factors. (those only make sense for the
82  // scaled points, for the unscaled ones they are 1.0.)
83  Scalar pcowLeverettFactor;
84  Scalar pcgoLeverettFactor;
85 
86  // Scaled relative permeabilities at residual displacing saturation
87  Scalar Krwr; // water
88  Scalar Krgr; // gas
89  Scalar Krorw; // oil in water-oil system
90  Scalar Krorg; // oil in gas-oil system
91 
92  // maximum relative permabilities
93  Scalar maxKrw; // maximum relative permability of water
94  Scalar maxKrow; // maximum relative permability of oil in the oil-water system
95  Scalar maxKrog; // maximum relative permability of oil in the gas-oil system
96  Scalar maxKrg; // maximum relative permability of gas
97 
98  bool operator==(const EclEpsScalingPointsInfo<Scalar>& data) const
99  {
100  return Swl == data.Swl &&
101  Sgl == data.Sgl &&
102  Swcr == data.Swcr &&
103  Sgcr == data.Sgcr &&
104  Sowcr == data.Sowcr &&
105  Sogcr == data.Sogcr &&
106  Swu == data.Swu &&
107  Sgu == data.Sgu &&
108  maxPcow == data.maxPcow &&
109  maxPcgo == data.maxPcgo &&
110  pcowLeverettFactor == data.pcowLeverettFactor &&
111  pcgoLeverettFactor == data.pcgoLeverettFactor &&
112  Krwr == data.Krwr &&
113  Krgr == data.Krgr &&
114  Krorw == data.Krorw &&
115  Krorg == data.Krorg &&
116  maxKrw == data.maxKrw &&
117  maxKrow == data.maxKrow &&
118  maxKrog == data.maxKrog &&
119  maxKrg == data.maxKrg;
120  }
121 
122  void print() const
123  {
124  std::cout << " Swl: " << Swl << '\n'
125  << " Sgl: " << Sgl << '\n'
126  << " Swcr: " << Swcr << '\n'
127  << " Sgcr: " << Sgcr << '\n'
128  << " Sowcr: " << Sowcr << '\n'
129  << " Sogcr: " << Sogcr << '\n'
130  << " Swu: " << Swu << '\n'
131  << " Sgu: " << Sgu << '\n'
132  << " maxPcow: " << maxPcow << '\n'
133  << " maxPcgo: " << maxPcgo << '\n'
134  << " pcowLeverettFactor: " << pcowLeverettFactor << '\n'
135  << " pcgoLeverettFactor: " << pcgoLeverettFactor << '\n'
136  << " Krwr: " << Krwr << '\n'
137  << " Krgr: " << Krgr << '\n'
138  << " Krorw: " << Krorw << '\n'
139  << " Krorg: " << Krorg << '\n'
140  << " maxKrw: " << maxKrw << '\n'
141  << " maxKrg: " << maxKrg << '\n'
142  << " maxKrow: " << maxKrow << '\n'
143  << " maxKrog: " << maxKrog << '\n';
144  }
145 
146 #if HAVE_ECL_INPUT
153  void extractUnscaled(const satfunc::RawTableEndPoints& rtep,
154  const satfunc::RawFunctionValues& rfunc,
155  const std::vector<double>::size_type satRegionIdx)
156  {
157  this->Swl = rtep.connate.water[satRegionIdx];
158  this->Sgl = rtep.connate.gas [satRegionIdx];
159 
160  this->Swcr = rtep.critical.water [satRegionIdx];
161  this->Sgcr = rtep.critical.gas [satRegionIdx];
162  this->Sowcr = rtep.critical.oil_in_water[satRegionIdx];
163  this->Sogcr = rtep.critical.oil_in_gas [satRegionIdx];
164 
165  this->Swu = rtep.maximum.water[satRegionIdx];
166  this->Sgu = rtep.maximum.gas [satRegionIdx];
167 
168  this->maxPcgo = rfunc.pc.g[satRegionIdx];
169  this->maxPcow = rfunc.pc.w[satRegionIdx];
170 
171  // there are no "unscaled" Leverett factors, so we just set them to 1.0
172  this->pcowLeverettFactor = 1.0;
173  this->pcgoLeverettFactor = 1.0;
174 
175  this->Krwr = rfunc.krw.r [satRegionIdx];
176  this->Krgr = rfunc.krg.r [satRegionIdx];
177  this->Krorw = rfunc.kro.rw[satRegionIdx];
178  this->Krorg = rfunc.kro.rg[satRegionIdx];
179 
180  this->maxKrw = rfunc.krw.max[satRegionIdx];
181  this->maxKrow = rfunc.kro.max[satRegionIdx];
182  this->maxKrog = rfunc.kro.max[satRegionIdx];
183  this->maxKrg = rfunc.krg.max[satRegionIdx];
184  }
185 
186  void update(Scalar& targetValue, const double * value_ptr) {
187  if (value_ptr)
188  targetValue = *value_ptr;
189  }
190 
196  void extractScaled(const EclipseState& eclState,
197  const EclEpsGridProperties& epsProperties,
198  unsigned activeIndex)
199  {
200  // overwrite the unscaled values with the values for the cell if it is
201  // explicitly specified by the corresponding keyword.
202  update(Swl, epsProperties.swl(activeIndex));
203  update(Sgl, epsProperties.sgl(activeIndex));
204  update(Swcr, epsProperties.swcr(activeIndex));
205  update(Sgcr, epsProperties.sgcr(activeIndex));
206 
207  update(Sowcr, epsProperties.sowcr(activeIndex));
208  update(Sogcr, epsProperties.sogcr(activeIndex));
209  update(Swu, epsProperties.swu(activeIndex));
210  update(Sgu, epsProperties.sgu(activeIndex));
211  update(maxPcow, epsProperties.pcw(activeIndex));
212  update(maxPcgo, epsProperties.pcg(activeIndex));
213 
214  update(this->Krwr, epsProperties.krwr(activeIndex));
215  update(this->Krgr, epsProperties.krgr(activeIndex));
216  update(this->Krorw, epsProperties.krorw(activeIndex));
217  update(this->Krorg, epsProperties.krorg(activeIndex));
218 
219  update(maxKrw, epsProperties.krw(activeIndex));
220  update(maxKrg, epsProperties.krg(activeIndex));
221  update(maxKrow, epsProperties.kro(activeIndex));
222  update(maxKrog, epsProperties.kro(activeIndex));
223 
224  // compute the Leverett capillary pressure scaling factors if applicable. note
225  // that this needs to be done using non-SI units to make it correspond to the
226  // documentation.
227  pcowLeverettFactor = 1.0;
228  pcgoLeverettFactor = 1.0;
229  if (eclState.getTableManager().useJFunc()) {
230  const auto& jfunc = eclState.getTableManager().getJFunc();
231  const auto& jfuncDir = jfunc.direction();
232 
233  Scalar perm;
234  if (jfuncDir == JFunc::Direction::X)
235  perm = epsProperties.permx(activeIndex);
236  else if (jfuncDir == JFunc::Direction::Y)
237  perm = epsProperties.permy(activeIndex);
238  else if (jfuncDir == JFunc::Direction::Z)
239  perm = epsProperties.permz(activeIndex);
240  else if (jfuncDir == JFunc::Direction::XY)
241  // TODO: verify that this really is the arithmetic mean. (the
242  // documentation just says that the "average" should be used, IMO the
243  // harmonic mean would be more appropriate because that's what's usually
244  // applied when calculating the fluxes.)
245  {
246  double permx = epsProperties.permx(activeIndex);
247  double permy = epsProperties.permy(activeIndex);
248  perm = arithmeticMean(permx, permy);
249  } else
250  throw std::runtime_error("Illegal direction indicator for the JFUNC "
251  "keyword ("+std::to_string(int(jfuncDir))+")");
252 
253  // convert permeability from m^2 to mD
254  perm *= 1.01325e15;
255 
256  Scalar poro = epsProperties.poro(activeIndex);
257  Scalar alpha = jfunc.alphaFactor();
258  Scalar beta = jfunc.betaFactor();
259 
260  // the part of the Leverett capillary pressure which does not depend on
261  // surface tension.
262  Scalar commonFactor = std::pow(poro, alpha)/std::pow(perm, beta);
263 
264  // multiply the documented constant by 10^5 because we want the pressures
265  // in [Pa], not in [bar]
266  const Scalar Uconst = 0.318316 * 1e5;
267 
268  // compute the oil-water Leverett factor.
269  const auto& jfuncFlag = jfunc.flag();
270  if (jfuncFlag == JFunc::Flag::WATER || jfuncFlag == JFunc::Flag::BOTH) {
271  // note that we use the surface tension in terms of [dyn/cm]
272  Scalar gamma =
273  jfunc.owSurfaceTension();
274  pcowLeverettFactor = commonFactor*gamma*Uconst;
275  }
276 
277  // compute the gas-oil Leverett factor.
278  if (jfuncFlag == JFunc::Flag::GAS || jfuncFlag == JFunc::Flag::BOTH) {
279  // note that we use the surface tension in terms of [dyn/cm]
280  Scalar gamma =
281  jfunc.goSurfaceTension();
282  pcgoLeverettFactor = commonFactor*gamma*Uconst;
283  }
284  }
285  }
286 #endif
287 
288 private:
289  void extractGridPropertyValue_(Scalar& targetValue,
290  const std::vector<double>* propData,
291  unsigned cartesianCellIdx)
292  {
293  if (!propData)
294  return;
295 
296  targetValue = (*propData)[cartesianCellIdx];
297  }
298 };
299 
306 template <class Scalar>
308 {
309 public:
314  const EclEpsConfig& config,
315  EclTwoPhaseSystemType epsSystemType)
316  {
317  if (epsSystemType == EclOilWaterSystem) {
318  // saturation scaling for capillary pressure
319  saturationPcPoints_[0] = epsInfo.Swl;
320  saturationPcPoints_[2] = saturationPcPoints_[1] = epsInfo.Swu;
321 
322  // krw saturation scaling endpoints
323  saturationKrwPoints_[0] = epsInfo.Swcr;
324  saturationKrwPoints_[1] = 1.0 - epsInfo.Sowcr - epsInfo.Sgl;
325  saturationKrwPoints_[2] = epsInfo.Swu;
326 
327  // krn saturation scaling endpoints (with the non-wetting phase being oil).
328  // because opm-material specifies non-wetting phase relperms in terms of the
329  // wetting phase saturations, the code here uses 1 minus the values specified
330  // by the Eclipse TD and the order of the scaling points is reversed
331  saturationKrnPoints_[2] = 1.0 - epsInfo.Sowcr;
332  saturationKrnPoints_[1] = epsInfo.Swcr + epsInfo.Sgl;
333  saturationKrnPoints_[0] = epsInfo.Swl + epsInfo.Sgl;
334 
335  if (config.enableLeverettScaling())
336  maxPcnwOrLeverettFactor_ = epsInfo.pcowLeverettFactor;
337  else
338  maxPcnwOrLeverettFactor_ = epsInfo.maxPcow;
339 
340  Krwr_ = epsInfo.Krwr;
341  Krnr_ = epsInfo.Krorw;
342 
343  maxKrw_ = epsInfo.maxKrw;
344  maxKrn_ = epsInfo.maxKrow;
345  }
346  else {
347  assert((epsSystemType == EclGasOilSystem) ||(epsSystemType == EclGasWaterSystem) );
348 
349  // saturation scaling for capillary pressure
350  saturationPcPoints_[0] = 1.0 - epsInfo.Swl - epsInfo.Sgu;
351  saturationPcPoints_[2] = saturationPcPoints_[1] = 1.0 - epsInfo.Swl - epsInfo.Sgl;
352 
353  // krw saturation scaling endpoints
354  saturationKrwPoints_[0] = epsInfo.Sogcr;
355  saturationKrwPoints_[1] = 1.0 - epsInfo.Sgcr - epsInfo.Swl;
356  saturationKrwPoints_[2] = 1.0 - epsInfo.Swl - epsInfo.Sgl;
357 
358  // krn saturation scaling endpoints (with the non-wetting phase being gas).
359  //
360  // As opm-material specifies non-wetting phase relative
361  // permeabilities in terms of the wetting phase saturations, the
362  // code here uses (1-SWL) minus the values specified by the
363  // ECLIPSE TD and the order of the scaling points is reversed.
364  saturationKrnPoints_[2] = 1.0 - epsInfo.Swl - epsInfo.Sgcr;
365  saturationKrnPoints_[1] = epsInfo.Sogcr;
366  saturationKrnPoints_[0] = 1.0 - epsInfo.Swl - epsInfo.Sgu;
367 
368  if (config.enableLeverettScaling())
369  maxPcnwOrLeverettFactor_ = epsInfo.pcgoLeverettFactor;
370  else
371  maxPcnwOrLeverettFactor_ = epsInfo.maxPcgo;
372 
373  Krwr_ = epsInfo.Krorg;
374  Krnr_ = epsInfo.Krgr;
375 
376  maxKrw_ = epsInfo.maxKrog;
377  maxKrn_ = epsInfo.maxKrg;
378  }
379  }
380 
384  void setSaturationPcPoint(unsigned pointIdx, Scalar value)
385  { saturationPcPoints_[pointIdx] = value; }
386 
390  const std::array<Scalar, 3>& saturationPcPoints() const
391  { return saturationPcPoints_; }
392 
396  void setSaturationKrwPoint(unsigned pointIdx, Scalar value)
397  { saturationKrwPoints_[pointIdx] = value; }
398 
402  const std::array<Scalar, 3>& saturationKrwPoints() const
403  { return saturationKrwPoints_; }
404 
408  void setSaturationKrnPoint(unsigned pointIdx, Scalar value)
409  { saturationKrnPoints_[pointIdx] = value; }
410 
414  const std::array<Scalar, 3>& saturationKrnPoints() const
415  { return saturationKrnPoints_; }
416 
420  void setMaxPcnw(Scalar value)
421  { maxPcnwOrLeverettFactor_ = value; }
422 
426  Scalar maxPcnw() const
427  { return maxPcnwOrLeverettFactor_; }
428 
432  void setLeverettFactor(Scalar value)
433  { maxPcnwOrLeverettFactor_ = value; }
434 
438  Scalar leverettFactor() const
439  { return maxPcnwOrLeverettFactor_; }
440 
445  void setKrwr(Scalar value)
446  { this->Krwr_ = value; }
447 
452  Scalar krwr() const
453  { return this->Krwr_; }
454 
458  void setMaxKrw(Scalar value)
459  { maxKrw_ = value; }
460 
464  Scalar maxKrw() const
465  { return maxKrw_; }
466 
471  void setKrnr(Scalar value)
472  { this->Krnr_ = value; }
473 
478  Scalar krnr() const
479  { return this->Krnr_; }
480 
484  void setMaxKrn(Scalar value)
485  { maxKrn_ = value; }
486 
490  Scalar maxKrn() const
491  { return maxKrn_; }
492 
493  void print() const
494  {
495  std::cout << " saturationKrnPoints_[0]: " << saturationKrnPoints_[0] << "\n"
496  << " saturationKrnPoints_[1]: " << saturationKrnPoints_[1] << "\n"
497  << " saturationKrnPoints_[2]: " << saturationKrnPoints_[2] << "\n";
498  }
499 
500 private:
501  // Points used for vertical scaling of capillary pressure
502  Scalar maxPcnwOrLeverettFactor_;
503 
504  // Maximum wetting phase relative permability value.
505  Scalar maxKrw_;
506 
507  // Scaled wetting phase relative permeability value at residual
508  // saturation of non-wetting phase.
509  Scalar Krwr_;
510 
511  // Maximum non-wetting phase relative permability value
512  Scalar maxKrn_;
513 
514  // Scaled non-wetting phase relative permeability value at residual
515  // saturation of wetting phase.
516  Scalar Krnr_;
517 
518  // The the points used for saturation ("x-axis") scaling of capillary pressure
519  std::array<Scalar, 3> saturationPcPoints_;
520 
521  // The the points used for saturation ("x-axis") scaling of wetting phase relative permeability
522  std::array<Scalar, 3> saturationKrwPoints_;
523 
524  // The the points used for saturation ("x-axis") scaling of non-wetting phase relative permeability
525  std::array<Scalar, 3> saturationKrnPoints_;
526 };
527 
528 } // namespace Opm
529 
530 #endif
Specifies the configuration used by the endpoint scaling code.
EclTwoPhaseSystemType
Specified which fluids are involved in a given twophase material law for endpoint scaling.
Definition: EclEpsConfig.hpp:49
Implements some common averages.
Scalar arithmeticMean(Scalar x, Scalar y)
Computes the arithmetic average of two values.
Definition: Means.hpp:45
Specifies the configuration used by the endpoint scaling code.
Definition: EclEpsConfig.hpp:63
bool enableLeverettScaling() const
Returns whether the Leverett capillary pressure scaling is enabled.
Definition: EclEpsConfig.hpp:183
Definition: EclEpsGridProperties.hpp:76
Represents the points on the X and Y axis to be scaled if endpoint scaling is used.
Definition: EclEpsScalingPoints.hpp:308
const std::array< Scalar, 3 > & saturationPcPoints() const
Returns the points used for capillary pressure saturation scaling.
Definition: EclEpsScalingPoints.hpp:390
void setMaxKrn(Scalar value)
Sets the maximum wetting phase relative permeability.
Definition: EclEpsScalingPoints.hpp:484
Scalar krnr() const
Returns non-wetting phase relative permeability at residual saturation of wetting phase.
Definition: EclEpsScalingPoints.hpp:478
Scalar maxKrn() const
Returns the maximum wetting phase relative permeability.
Definition: EclEpsScalingPoints.hpp:490
void setKrnr(Scalar value)
Set non-wetting phase relative permeability at residual saturation of wetting phase.
Definition: EclEpsScalingPoints.hpp:471
const std::array< Scalar, 3 > & saturationKrnPoints() const
Returns the points used for non-wetting phase relperm saturation scaling.
Definition: EclEpsScalingPoints.hpp:414
void setMaxKrw(Scalar value)
Sets the maximum wetting phase relative permeability.
Definition: EclEpsScalingPoints.hpp:458
Scalar maxKrw() const
Returns the maximum wetting phase relative permeability.
Definition: EclEpsScalingPoints.hpp:464
void setSaturationKrwPoint(unsigned pointIdx, Scalar value)
Sets an saturation value for wetting-phase relperm saturation scaling.
Definition: EclEpsScalingPoints.hpp:396
void setMaxPcnw(Scalar value)
Sets the maximum capillary pressure.
Definition: EclEpsScalingPoints.hpp:420
Scalar krwr() const
Returns wetting-phase relative permeability at residual saturation of non-wetting phase.
Definition: EclEpsScalingPoints.hpp:452
Scalar maxPcnw() const
Returns the maximum capillary pressure.
Definition: EclEpsScalingPoints.hpp:426
void setKrwr(Scalar value)
Set wetting-phase relative permeability at residual saturation of non-wetting phase.
Definition: EclEpsScalingPoints.hpp:445
void setSaturationPcPoint(unsigned pointIdx, Scalar value)
Sets an saturation value for capillary pressure saturation scaling.
Definition: EclEpsScalingPoints.hpp:384
void setLeverettFactor(Scalar value)
Sets the Leverett scaling factor for capillary pressure.
Definition: EclEpsScalingPoints.hpp:432
const std::array< Scalar, 3 > & saturationKrwPoints() const
Returns the points used for wetting phase relperm saturation scaling.
Definition: EclEpsScalingPoints.hpp:402
void init(const EclEpsScalingPointsInfo< Scalar > &epsInfo, const EclEpsConfig &config, EclTwoPhaseSystemType epsSystemType)
Assigns the scaling points which actually ought to be used.
Definition: EclEpsScalingPoints.hpp:313
void setSaturationKrnPoint(unsigned pointIdx, Scalar value)
Sets an saturation value for non-wetting phase relperm saturation scaling.
Definition: EclEpsScalingPoints.hpp:408
Scalar leverettFactor() const
Returns the Leverett scaling factor for capillary pressure.
Definition: EclEpsScalingPoints.hpp:438
This structure represents all values which can be possibly used as scaling points in the endpoint sca...
Definition: EclEpsScalingPoints.hpp:62