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