My Project
EclMaterialLawManager.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 #if ! HAVE_ECL_INPUT
28 #error "Eclipse input support in opm-common is required to use the ECL material manager!"
29 #endif
30 
31 #ifndef OPM_ECL_MATERIAL_LAW_MANAGER_HPP
32 #define OPM_ECL_MATERIAL_LAW_MANAGER_HPP
33 
44 
45 #if HAVE_OPM_COMMON
46 #include <opm/common/OpmLog/OpmLog.hpp>
47 #endif
48 
49 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
50 #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
51 #include <opm/input/eclipse/EclipseState/Tables/TableColumn.hpp>
52 
53 #include <algorithm>
54 #include <cassert>
55 #include <memory>
56 #include <stdexcept>
57 #include <vector>
58 
59 namespace Opm {
60 
67 template <class TraitsT>
69 {
70 private:
71  typedef TraitsT Traits;
72  typedef typename Traits::Scalar Scalar;
73  enum { waterPhaseIdx = Traits::wettingPhaseIdx };
74  enum { oilPhaseIdx = Traits::nonWettingPhaseIdx };
75  enum { gasPhaseIdx = Traits::gasPhaseIdx };
76  enum { numPhases = Traits::numPhases };
77 
81 
82  // the two-phase material law which is defined on effective (unscaled) saturations
86 
87  typedef typename GasOilEffectiveTwoPhaseLaw::Params GasOilEffectiveTwoPhaseParams;
88  typedef typename OilWaterEffectiveTwoPhaseLaw::Params OilWaterEffectiveTwoPhaseParams;
89  typedef typename GasWaterEffectiveTwoPhaseLaw::Params GasWaterEffectiveTwoPhaseParams;
90 
91  // the two-phase material law which is defined on absolute (scaled) saturations
95  typedef typename GasOilEpsTwoPhaseLaw::Params GasOilEpsTwoPhaseParams;
96  typedef typename OilWaterEpsTwoPhaseLaw::Params OilWaterEpsTwoPhaseParams;
97  typedef typename GasWaterEpsTwoPhaseLaw::Params GasWaterEpsTwoPhaseParams;
98 
99  // the scaled two-phase material laws with hystersis
103  typedef typename GasOilTwoPhaseLaw::Params GasOilTwoPhaseHystParams;
104  typedef typename OilWaterTwoPhaseLaw::Params OilWaterTwoPhaseHystParams;
105  typedef typename GasWaterTwoPhaseLaw::Params GasWaterTwoPhaseHystParams;
106 
107 public:
108  // the three-phase material law used by the simulation
110  typedef typename MaterialLaw::Params MaterialLawParams;
111 
112 private:
113  // internal typedefs
114  typedef std::vector<std::shared_ptr<GasOilEffectiveTwoPhaseParams> > GasOilEffectiveParamVector;
115  typedef std::vector<std::shared_ptr<OilWaterEffectiveTwoPhaseParams> > OilWaterEffectiveParamVector;
116  typedef std::vector<std::shared_ptr<GasWaterEffectiveTwoPhaseParams> > GasWaterEffectiveParamVector;
117 
118  typedef std::vector<std::shared_ptr<EclEpsScalingPoints<Scalar> > > GasOilScalingPointsVector;
119  typedef std::vector<std::shared_ptr<EclEpsScalingPoints<Scalar> > > OilWaterScalingPointsVector;
120  typedef std::vector<std::shared_ptr<EclEpsScalingPoints<Scalar> > > GasWaterScalingPointsVector;
121  typedef std::vector<EclEpsScalingPointsInfo<Scalar>> OilWaterScalingInfoVector;
122  typedef std::vector<std::shared_ptr<GasOilTwoPhaseHystParams> > GasOilParamVector;
123  typedef std::vector<std::shared_ptr<OilWaterTwoPhaseHystParams> > OilWaterParamVector;
124  typedef std::vector<std::shared_ptr<GasWaterTwoPhaseHystParams> > GasWaterParamVector;
125  typedef std::vector<std::shared_ptr<MaterialLawParams> > MaterialLawParamsVector;
126 
127 public:
129  {}
130 
131  void initFromState(const EclipseState& eclState)
132  {
133  // get the number of saturation regions and the number of cells in the deck
134  const auto& runspec = eclState.runspec();
135  const size_t numSatRegions = runspec.tabdims().getNumSatTables();
136 
137  const auto& ph = runspec.phases();
138  this->hasGas = ph.active(Phase::GAS);
139  this->hasOil = ph.active(Phase::OIL);
140  this->hasWater = ph.active(Phase::WATER);
141 
142  readGlobalEpsOptions_(eclState);
143  readGlobalHysteresisOptions_(eclState);
144  readGlobalThreePhaseOptions_(runspec);
145 
146  // Read the end point scaling configuration (once per run).
147  gasOilConfig = std::make_shared<EclEpsConfig>();
148  oilWaterConfig = std::make_shared<EclEpsConfig>();
149  gasWaterConfig = std::make_shared<EclEpsConfig>();
150  gasOilConfig->initFromState(eclState, EclGasOilSystem);
151  oilWaterConfig->initFromState(eclState, EclOilWaterSystem);
152  gasWaterConfig->initFromState(eclState, EclGasWaterSystem);
153 
154 
155  const auto& tables = eclState.getTableManager();
156 
157  {
158  const auto& stone1exTables = tables.getStone1exTable();
159 
160  if (! stone1exTables.empty()) {
161  stoneEtas.clear();
162  stoneEtas.reserve(numSatRegions);
163 
164  for (const auto& table : stone1exTables) {
165  stoneEtas.push_back(table.eta);
166  }
167  }
168  }
169 
170  this->unscaledEpsInfo_.resize(numSatRegions);
171 
172  if (this->hasGas + this->hasOil + this->hasWater == 1) {
173  // Single-phase simulation. Special case. Nothing to do here.
174  return;
175  }
176 
177  // Multiphase simulation. Common case.
178  const auto tolcrit = runspec.saturationFunctionControls()
179  .minimumRelpermMobilityThreshold();
180 
181  const auto rtep = satfunc::getRawTableEndpoints(tables, ph, tolcrit);
182  const auto rfunc = satfunc::getRawFunctionValues(tables, ph, rtep);
183 
184  for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) {
185  this->unscaledEpsInfo_[satRegionIdx]
186  .extractUnscaled(rtep, rfunc, satRegionIdx);
187  }
188  }
189 
190  void initParamsForElements(const EclipseState& eclState, size_t numCompressedElems)
191  {
192  // get the number of saturation regions
193  const size_t numSatRegions = eclState.runspec().tabdims().getNumSatTables();
194 
195  // setup the saturation region specific parameters
196  gasOilUnscaledPointsVector_.resize(numSatRegions);
197  oilWaterUnscaledPointsVector_.resize(numSatRegions);
198  gasWaterUnscaledPointsVector_.resize(numSatRegions);
199 
200  gasOilEffectiveParamVector_.resize(numSatRegions);
201  oilWaterEffectiveParamVector_.resize(numSatRegions);
202  gasWaterEffectiveParamVector_.resize(numSatRegions);
203  for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) {
204  // unscaled points for end-point scaling
205  readGasOilUnscaledPoints_(gasOilUnscaledPointsVector_, gasOilConfig, eclState, satRegionIdx);
206  readOilWaterUnscaledPoints_(oilWaterUnscaledPointsVector_, oilWaterConfig, eclState, satRegionIdx);
207  readGasWaterUnscaledPoints_(gasWaterUnscaledPointsVector_, gasWaterConfig, eclState, satRegionIdx);
208 
209  // the parameters for the effective two-phase matererial laws
210  readGasOilEffectiveParameters_(gasOilEffectiveParamVector_, eclState, satRegionIdx);
211  readOilWaterEffectiveParameters_(oilWaterEffectiveParamVector_, eclState, satRegionIdx);
212  readGasWaterEffectiveParameters_(gasWaterEffectiveParamVector_, eclState, satRegionIdx);
213  }
214 
215  // copy the SATNUM grid property. in some cases this is not necessary, but it
216  // should not require much memory anyway...
217  satnumRegionArray_.resize(numCompressedElems);
218  if (eclState.fieldProps().has_int("SATNUM")) {
219  const auto& satnumRawData = eclState.fieldProps().get_int("SATNUM");
220  for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
221  satnumRegionArray_[elemIdx] = satnumRawData[elemIdx] - 1;
222  }
223  }
224  else
225  std::fill(satnumRegionArray_.begin(), satnumRegionArray_.end(), 0);
226 
227  // create the information for the imbibition region (IMBNUM). By default this is
228  // the same as the saturation region (SATNUM)
229  imbnumRegionArray_ = satnumRegionArray_;
230  if (eclState.fieldProps().has_int("IMBNUM")) {
231  const auto& imbnumRawData = eclState.fieldProps().get_int("IMBNUM");
232  for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
233  imbnumRegionArray_[elemIdx] = imbnumRawData[elemIdx] - 1;
234  }
235  }
236 
237  assert(numCompressedElems == satnumRegionArray_.size());
238  assert(!enableHysteresis() || numCompressedElems == imbnumRegionArray_.size());
239 
240  // read the scaled end point scaling parameters which are specific for each
241  // element
242  oilWaterScaledEpsInfoDrainage_.resize(numCompressedElems);
243 
244  std::unique_ptr<EclEpsGridProperties> epsImbGridProperties;
245 
246  if (enableHysteresis()) {
247  epsImbGridProperties = std::make_unique<EclEpsGridProperties>(eclState, true);
248  }
249 
250  EclEpsGridProperties epsGridProperties(eclState, false);
251  materialLawParams_.resize(numCompressedElems);
252 
253  for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
254  unsigned satRegionIdx = static_cast<unsigned>(satnumRegionArray_[elemIdx]);
255  auto gasOilParams = std::make_shared<GasOilTwoPhaseHystParams>();
256  auto oilWaterParams = std::make_shared<OilWaterTwoPhaseHystParams>();
257  auto gasWaterParams = std::make_shared<GasWaterTwoPhaseHystParams>();
258  gasOilParams->setConfig(hysteresisConfig_);
259  oilWaterParams->setConfig(hysteresisConfig_);
260  gasWaterParams->setConfig(hysteresisConfig_);
261 
262  auto [gasOilScaledInfo, gasOilScaledPoint] =
263  readScaledPoints_(*gasOilConfig,
264  eclState,
265  epsGridProperties,
266  elemIdx,
267  EclGasOilSystem);
268 
269  auto [owinfo, oilWaterScaledEpsPointDrainage] =
270  readScaledPoints_(*oilWaterConfig,
271  eclState,
272  epsGridProperties,
273  elemIdx,
274  EclOilWaterSystem);
275  oilWaterScaledEpsInfoDrainage_[elemIdx] = owinfo;
276 
277  auto [gasWaterScaledInfo, gasWaterScaledPoint] =
278  readScaledPoints_(*gasWaterConfig,
279  eclState,
280  epsGridProperties,
281  elemIdx,
282  EclGasWaterSystem);
283 
284  if (hasGas && hasOil) {
285  GasOilEpsTwoPhaseParams gasOilDrainParams;
286  gasOilDrainParams.setConfig(gasOilConfig);
287  gasOilDrainParams.setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
288  gasOilDrainParams.setScaledPoints(gasOilScaledPoint);
289  gasOilDrainParams.setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
290  gasOilDrainParams.finalize();
291 
292  gasOilParams->setDrainageParams(gasOilDrainParams);
293  }
294 
295  if (hasOil && hasWater) {
296  OilWaterEpsTwoPhaseParams oilWaterDrainParams;
297  oilWaterDrainParams.setConfig(oilWaterConfig);
298  oilWaterDrainParams.setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
299  oilWaterDrainParams.setScaledPoints(oilWaterScaledEpsPointDrainage);
300  oilWaterDrainParams.setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
301  oilWaterDrainParams.finalize();
302 
303  oilWaterParams->setDrainageParams(oilWaterDrainParams);
304  }
305 
306  if (hasGas && hasWater && !hasOil) {
307  GasWaterEpsTwoPhaseParams gasWaterDrainParams;
308  gasWaterDrainParams.setConfig(gasWaterConfig);
309  gasWaterDrainParams.setUnscaledPoints(gasWaterUnscaledPointsVector_[satRegionIdx]);
310  gasWaterDrainParams.setScaledPoints(gasWaterScaledPoint);
311  gasWaterDrainParams.setEffectiveLawParams(gasWaterEffectiveParamVector_[satRegionIdx]);
312  gasWaterDrainParams.finalize();
313 
314  gasWaterParams->setDrainageParams(gasWaterDrainParams);
315  }
316 
317  if (enableHysteresis()) {
318  auto [gasOilScaledImbInfo, gasOilScaledImbPoint] =
319  readScaledPoints_(*gasOilConfig,
320  eclState,
321  *epsImbGridProperties,
322  elemIdx,
323  EclGasOilSystem);
324 
325  auto [oilWaterScaledImbInfo, oilWaterScaledImbPoint] =
326  readScaledPoints_(*oilWaterConfig,
327  eclState,
328  *epsImbGridProperties,
329  elemIdx,
330  EclOilWaterSystem);
331 
332  auto [gasWaterScaledImbInfo, gasWaterScaledImbPoint] =
333  readScaledPoints_(*gasWaterConfig,
334  eclState,
335  *epsImbGridProperties,
336  elemIdx,
337  EclGasWaterSystem);
338 
339  unsigned imbRegionIdx = imbnumRegionArray_[elemIdx];
340  if (hasGas && hasOil) {
341  GasOilEpsTwoPhaseParams gasOilImbParamsHyst;
342  gasOilImbParamsHyst.setConfig(gasOilConfig);
343  gasOilImbParamsHyst.setUnscaledPoints(gasOilUnscaledPointsVector_[imbRegionIdx]);
344  gasOilImbParamsHyst.setScaledPoints(gasOilScaledImbPoint);
345  gasOilImbParamsHyst.setEffectiveLawParams(gasOilEffectiveParamVector_[imbRegionIdx]);
346  gasOilImbParamsHyst.finalize();
347 
348  gasOilParams->setImbibitionParams(gasOilImbParamsHyst,
349  gasOilScaledImbInfo,
350  EclGasOilSystem);
351  }
352 
353  if (hasOil && hasWater) {
354  OilWaterEpsTwoPhaseParams oilWaterImbParamsHyst;
355  oilWaterImbParamsHyst.setConfig(oilWaterConfig);
356  oilWaterImbParamsHyst.setUnscaledPoints(oilWaterUnscaledPointsVector_[imbRegionIdx]);
357  oilWaterImbParamsHyst.setScaledPoints(oilWaterScaledImbPoint);
358  oilWaterImbParamsHyst.setEffectiveLawParams(oilWaterEffectiveParamVector_[imbRegionIdx]);
359  oilWaterImbParamsHyst.finalize();
360 
361  oilWaterParams->setImbibitionParams(oilWaterImbParamsHyst,
362  gasOilScaledImbInfo,
363  EclOilWaterSystem);
364  }
365 
366  if (hasGas && hasWater && !hasOil) {
367  GasWaterEpsTwoPhaseParams gasWaterImbParamsHyst;
368  gasWaterImbParamsHyst.setConfig(gasWaterConfig);
369  gasWaterImbParamsHyst.setUnscaledPoints(gasWaterUnscaledPointsVector_[imbRegionIdx]);
370  gasWaterImbParamsHyst.setScaledPoints(gasWaterScaledImbPoint);
371  gasWaterImbParamsHyst.setEffectiveLawParams(gasWaterEffectiveParamVector_[imbRegionIdx]);
372  gasWaterImbParamsHyst.finalize();
373 
374  gasWaterParams->setImbibitionParams(gasWaterImbParamsHyst,
375  gasWaterScaledImbInfo,
376  EclGasWaterSystem);
377  }
378  }
379 
380  if (hasGas && hasOil)
381  gasOilParams->finalize();
382 
383  if (hasOil && hasWater)
384  oilWaterParams->finalize();
385 
386  if (hasGas && hasWater && !hasOil)
387  gasWaterParams->finalize();
388 
389  initThreePhaseParams_(eclState,
390  materialLawParams_[elemIdx],
391  satRegionIdx,
392  oilWaterScaledEpsInfoDrainage_[elemIdx],
393  oilWaterParams,
394  gasOilParams,
395  gasWaterParams);
396 
397  materialLawParams_[elemIdx].finalize();
398  }
399  }
400 
401 
410  Scalar applySwatinit(unsigned elemIdx,
411  Scalar pcow,
412  Scalar Sw)
413  {
414  auto& elemScaledEpsInfo = oilWaterScaledEpsInfoDrainage_[elemIdx];
415 
416  // TODO: Mixed wettability systems - see ecl kw OPTIONS switch 74
417 
418  if (pcow < 0.0)
419  Sw = elemScaledEpsInfo.Swu;
420  else {
421 
422  if (Sw <= elemScaledEpsInfo.Swl)
423  Sw = elemScaledEpsInfo.Swl;
424 
425  // specify a fluid state which only stores the saturations
426  typedef SimpleModularFluidState<Scalar,
427  numPhases,
428  /*numComponents=*/0,
429  /*FluidSystem=*/void, /* -> don't care */
430  /*storePressure=*/false,
431  /*storeTemperature=*/false,
432  /*storeComposition=*/false,
433  /*storeFugacity=*/false,
434  /*storeSaturation=*/true,
435  /*storeDensity=*/false,
436  /*storeViscosity=*/false,
437  /*storeEnthalpy=*/false> FluidState;
438  FluidState fs;
439  fs.setSaturation(waterPhaseIdx, Sw);
440  fs.setSaturation(gasPhaseIdx, 0);
441  fs.setSaturation(oilPhaseIdx, 0);
442  Scalar pc[numPhases] = { 0 };
443  MaterialLaw::capillaryPressures(pc, materialLawParams(elemIdx), fs);
444 
445  Scalar pcowAtSw = pc[oilPhaseIdx] - pc[waterPhaseIdx];
446  const Scalar pcowAtSwThreshold = 1.0; //Pascal
447  // avoid divison by very small number
448  if (std::abs(pcowAtSw) > pcowAtSwThreshold) {
449  elemScaledEpsInfo.maxPcow *= pcow/pcowAtSw;
450  auto& elemEclEpsScalingPoints = oilWaterScaledEpsPointsDrainage(elemIdx);
451  elemEclEpsScalingPoints.init(elemScaledEpsInfo, *oilWaterEclEpsConfig_, EclOilWaterSystem);
452  }
453  }
454 
455  return Sw;
456  }
457 
458  bool enableEndPointScaling() const
459  { return enableEndPointScaling_; }
460 
461  bool enableHysteresis() const
462  { return hysteresisConfig_->enableHysteresis(); }
463 
464  MaterialLawParams& materialLawParams(unsigned elemIdx)
465  {
466  assert(elemIdx < materialLawParams_.size());
467  return materialLawParams_[elemIdx];
468  }
469 
470  const MaterialLawParams& materialLawParams(unsigned elemIdx) const
471  {
472  assert(elemIdx < materialLawParams_.size());
473  return materialLawParams_[elemIdx];
474  }
475 
484  const MaterialLawParams& connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const
485  {
486  MaterialLawParams& mlp = const_cast<MaterialLawParams&>(materialLawParams_[elemIdx]);
487 
488 #if HAVE_OPM_COMMON
489  if (enableHysteresis())
490  OpmLog::warning("Warning: Using non-default satnum regions for connection is not tested in combination with hysteresis");
491 #endif
492  // Currently we don't support COMPIMP. I.e. use the same table lookup for the hysteresis curves.
493  // unsigned impRegionIdx = satRegionIdx;
494 
495  // change the sat table it points to.
496  switch (mlp.approach()) {
497  case EclMultiplexerApproach::EclStone1Approach: {
498  auto& realParams = mlp.template getRealParams<EclMultiplexerApproach::EclStone1Approach>();
499 
500  realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
501  realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
502  realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
503  realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
504 // if (enableHysteresis()) {
505 // realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]);
506 // realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]);
507 // realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]);
508 // realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]);
509 // }
510  }
511  break;
512 
513  case EclMultiplexerApproach::EclStone2Approach: {
514  auto& realParams = mlp.template getRealParams<EclMultiplexerApproach::EclStone2Approach>();
515  realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
516  realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
517  realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
518  realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
519 // if (enableHysteresis()) {
520 // realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]);
521 // realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]);
522 // realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]);
523 // realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]);
524 // }
525  }
526  break;
527 
528  case EclMultiplexerApproach::EclDefaultApproach: {
529  auto& realParams = mlp.template getRealParams<EclMultiplexerApproach::EclDefaultApproach>();
530  realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
531  realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
532  realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
533  realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
534 // if (enableHysteresis()) {
535 // realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]);
536 // realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]);
537 // realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]);
538 // realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]);
539 // }
540  }
541  break;
542 
543  case EclMultiplexerApproach::EclTwoPhaseApproach: {
544  auto& realParams = mlp.template getRealParams<EclMultiplexerApproach::EclTwoPhaseApproach>();
545  realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
546  realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
547  realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
548  realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
549 // if (enableHysteresis()) {
550 // realParams.oilWaterParams().imbibitionParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[impRegionIdx]);
551 // realParams.oilWaterParams().imbibitionParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[impRegionIdx]);
552 // realParams.gasOilParams().imbibitionParams().setUnscaledPoints(gasOilUnscaledPointsVector_[impRegionIdx]);
553 // realParams.gasOilParams().imbibitionParams().setEffectiveLawParams(gasOilEffectiveParamVector_[impRegionIdx]);
554 // }
555  }
556  break;
557 
558  default:
559  throw std::logic_error("Enum value for material approach unknown!");
560  }
561 
562  return mlp;
563  }
564 
565  int satnumRegionIdx(unsigned elemIdx) const
566  { return satnumRegionArray_[elemIdx]; }
567 
568  int imbnumRegionIdx(unsigned elemIdx) const
569  { return imbnumRegionArray_[elemIdx]; }
570 
571  std::shared_ptr<MaterialLawParams>& materialLawParamsPointerReferenceHack(unsigned elemIdx)
572  {
573  assert(0 <= elemIdx && elemIdx < materialLawParams_.size());
574  return materialLawParams_[elemIdx];
575  }
576 
577  template <class FluidState>
578  void updateHysteresis(const FluidState& fluidState, unsigned elemIdx)
579  {
580  if (!enableHysteresis())
581  return;
582 
583  MaterialLaw::updateHysteresis(materialLawParams_[elemIdx], fluidState);
584  }
585 
586  void oilWaterHysteresisParams(Scalar& pcSwMdc,
587  Scalar& krnSwMdc,
588  unsigned elemIdx) const
589  {
590  if (!enableHysteresis())
591  throw std::runtime_error("Cannot get hysteresis parameters if hysteresis not enabled.");
592 
593  const auto& params = materialLawParams(elemIdx);
594  MaterialLaw::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, params);
595  }
596 
597  void setOilWaterHysteresisParams(const Scalar& pcSwMdc,
598  const Scalar& krnSwMdc,
599  unsigned elemIdx)
600  {
601  if (!enableHysteresis())
602  throw std::runtime_error("Cannot set hysteresis parameters if hysteresis not enabled.");
603 
604  auto& params = materialLawParams(elemIdx);
605  MaterialLaw::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, params);
606  }
607 
608  void gasOilHysteresisParams(Scalar& pcSwMdc,
609  Scalar& krnSwMdc,
610  unsigned elemIdx) const
611  {
612  if (!enableHysteresis())
613  throw std::runtime_error("Cannot get hysteresis parameters if hysteresis not enabled.");
614 
615  const auto& params = materialLawParams(elemIdx);
616  MaterialLaw::gasOilHysteresisParams(pcSwMdc, krnSwMdc, params);
617  }
618 
619  void setGasOilHysteresisParams(const Scalar& pcSwMdc,
620  const Scalar& krnSwMdc,
621  unsigned elemIdx)
622  {
623  if (!enableHysteresis())
624  throw std::runtime_error("Cannot set hysteresis parameters if hysteresis not enabled.");
625 
626  auto& params = materialLawParams(elemIdx);
627  MaterialLaw::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, params);
628  }
629 
630  EclEpsScalingPoints<Scalar>& oilWaterScaledEpsPointsDrainage(unsigned elemIdx)
631  {
632  auto& materialParams = materialLawParams_[elemIdx];
633  switch (materialParams.approach()) {
634  case EclMultiplexerApproach::EclStone1Approach: {
635  auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclStone1Approach>();
636  return realParams.oilWaterParams().drainageParams().scaledPoints();
637  }
638 
639  case EclMultiplexerApproach::EclStone2Approach: {
640  auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclStone2Approach>();
641  return realParams.oilWaterParams().drainageParams().scaledPoints();
642  }
643 
644  case EclMultiplexerApproach::EclDefaultApproach: {
645  auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclDefaultApproach>();
646  return realParams.oilWaterParams().drainageParams().scaledPoints();
647  }
648 
649  case EclMultiplexerApproach::EclTwoPhaseApproach: {
650  auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclTwoPhaseApproach>();
651  return realParams.oilWaterParams().drainageParams().scaledPoints();
652  }
653  default:
654  throw std::logic_error("Enum value for material approach unknown!");
655  }
656  }
657 
658  const EclEpsScalingPointsInfo<Scalar>& oilWaterScaledEpsInfoDrainage(size_t elemIdx) const
659  { return oilWaterScaledEpsInfoDrainage_[elemIdx]; }
660 
661 private:
662  void readGlobalEpsOptions_(const EclipseState& eclState)
663  {
664  oilWaterEclEpsConfig_ = std::make_shared<EclEpsConfig>();
665  oilWaterEclEpsConfig_->initFromState(eclState, EclOilWaterSystem);
666 
667  enableEndPointScaling_ = eclState.getTableManager().hasTables("ENKRVD");
668  }
669 
670  void readGlobalHysteresisOptions_(const EclipseState& state)
671  {
672  hysteresisConfig_ = std::make_shared<EclHysteresisConfig>();
673  hysteresisConfig_->initFromState(state.runspec());
674  }
675 
676  void readGlobalThreePhaseOptions_(const Runspec& runspec)
677  {
678  bool gasEnabled = runspec.phases().active(Phase::GAS);
679  bool oilEnabled = runspec.phases().active(Phase::OIL);
680  bool waterEnabled = runspec.phases().active(Phase::WATER);
681 
682  int numEnabled =
683  (gasEnabled?1:0)
684  + (oilEnabled?1:0)
685  + (waterEnabled?1:0);
686 
687  if (numEnabled == 0) {
688  throw std::runtime_error("At least one fluid phase must be enabled. (Is: "+std::to_string(numEnabled)+")");
689  } else if (numEnabled == 1) {
690  threePhaseApproach_ = EclMultiplexerApproach::EclOnePhaseApproach;
691  } else if ( numEnabled == 2) {
692  threePhaseApproach_ = EclMultiplexerApproach::EclTwoPhaseApproach;
693  if (!gasEnabled)
694  twoPhaseApproach_ = EclTwoPhaseApproach::EclTwoPhaseOilWater;
695  else if (!oilEnabled)
696  twoPhaseApproach_ = EclTwoPhaseApproach::EclTwoPhaseGasWater;
697  else if (!waterEnabled)
698  twoPhaseApproach_ = EclTwoPhaseApproach::EclTwoPhaseGasOil;
699  }
700  else {
701  assert(numEnabled == 3);
702 
703  threePhaseApproach_ = EclMultiplexerApproach::EclDefaultApproach;
704  const auto& satctrls = runspec.saturationFunctionControls();
705  if (satctrls.krModel() == SatFuncControls::ThreePhaseOilKrModel::Stone2)
706  threePhaseApproach_ = EclMultiplexerApproach::EclStone2Approach;
707  else if (satctrls.krModel() == SatFuncControls::ThreePhaseOilKrModel::Stone1)
708  threePhaseApproach_ = EclMultiplexerApproach::EclStone1Approach;
709  }
710  }
711 
712  template <class Container>
713  void readGasOilEffectiveParameters_(Container& dest,
714  const EclipseState& eclState,
715  unsigned satRegionIdx)
716  {
717  if (!hasGas || !hasOil)
718  // we don't read anything if either the gas or the oil phase is not active
719  return;
720 
721  dest[satRegionIdx] = std::make_shared<GasOilEffectiveTwoPhaseParams>();
722 
723  auto& effParams = *dest[satRegionIdx];
724 
725  // the situation for the gas phase is complicated that all saturations are
726  // shifted by the connate water saturation.
727  const Scalar Swco = unscaledEpsInfo_[satRegionIdx].Swl;
728  const auto tolcrit = eclState.runspec().saturationFunctionControls()
729  .minimumRelpermMobilityThreshold();
730 
731  const auto& tableManager = eclState.getTableManager();
732 
733  switch (eclState.runspec().saturationFunctionControls().family()) {
734  case SatFuncControls::KeywordFamily::Family_I:
735  {
736  const TableContainer& sgofTables = tableManager.getSgofTables();
737  const TableContainer& slgofTables = tableManager.getSlgofTables();
738  if (!sgofTables.empty())
739  readGasOilEffectiveParametersSgof_(effParams, Swco, tolcrit,
740  sgofTables.getTable<SgofTable>(satRegionIdx));
741  else if (!slgofTables.empty())
742  readGasOilEffectiveParametersSlgof_(effParams, Swco, tolcrit,
743  slgofTables.getTable<SlgofTable>(satRegionIdx));
744  break;
745  }
746 
747  case SatFuncControls::KeywordFamily::Family_II:
748  {
749  const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable<SgfnTable>( satRegionIdx );
750  if (!hasWater) {
751  // oil and gas case
752  const Sof2Table& sof2Table = tableManager.getSof2Tables().getTable<Sof2Table>( satRegionIdx );
753  readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof2Table, sgfnTable);
754  }
755  else {
756  const Sof3Table& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>( satRegionIdx );
757  readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof3Table, sgfnTable);
758  }
759  break;
760  }
761 
762  case SatFuncControls::KeywordFamily::Undefined:
763  throw std::domain_error("No valid saturation keyword family specified");
764  }
765  }
766 
767  void readGasOilEffectiveParametersSgof_(GasOilEffectiveTwoPhaseParams& effParams,
768  const Scalar Swco,
769  const double tolcrit,
770  const SgofTable& sgofTable)
771  {
772  // convert the saturations of the SGOF keyword from gas to oil saturations
773  std::vector<double> SoSamples(sgofTable.numRows());
774  for (size_t sampleIdx = 0; sampleIdx < sgofTable.numRows(); ++ sampleIdx) {
775  SoSamples[sampleIdx] = (1.0 - Swco) - sgofTable.get("SG", sampleIdx);
776  }
777 
778  effParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn("KROG")));
779  effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn("KRG")));
780  effParams.setPcnwSamples(SoSamples, sgofTable.getColumn("PCOG").vectorCopy());
781  effParams.finalize();
782  }
783 
784  void readGasOilEffectiveParametersSlgof_(GasOilEffectiveTwoPhaseParams& effParams,
785  const Scalar Swco,
786  const double tolcrit,
787  const SlgofTable& slgofTable)
788  {
789  // convert the saturations of the SLGOF keyword from "liquid" to oil saturations
790  std::vector<double> SoSamples(slgofTable.numRows());
791  for (size_t sampleIdx = 0; sampleIdx < slgofTable.numRows(); ++ sampleIdx) {
792  SoSamples[sampleIdx] = slgofTable.get("SL", sampleIdx) - Swco;
793  }
794 
795  effParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn("KROG")));
796  effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn("KRG")));
797  effParams.setPcnwSamples(SoSamples, slgofTable.getColumn("PCOG").vectorCopy());
798  effParams.finalize();
799  }
800 
801  void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
802  const Scalar Swco,
803  const double tolcrit,
804  const Sof3Table& sof3Table,
805  const SgfnTable& sgfnTable)
806  {
807  // convert the saturations of the SGFN keyword from gas to oil saturations
808  std::vector<double> SoSamples(sgfnTable.numRows());
809  std::vector<double> SoColumn = sof3Table.getColumn("SO").vectorCopy();
810  for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) {
811  SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx);
812  }
813 
814  effParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROG")));
815  effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG")));
816  effParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy());
817  effParams.finalize();
818  }
819 
820  void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
821  const Scalar Swco,
822  const double tolcrit,
823  const Sof2Table& sof2Table,
824  const SgfnTable& sgfnTable)
825  {
826  // convert the saturations of the SGFN keyword from gas to oil saturations
827  std::vector<double> SoSamples(sgfnTable.numRows());
828  std::vector<double> SoColumn = sof2Table.getColumn("SO").vectorCopy();
829  for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) {
830  SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get("SG", sampleIdx);
831  }
832 
833  effParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof2Table.getColumn("KRO")));
834  effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG")));
835  effParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy());
836  effParams.finalize();
837  }
838 
839  template <class Container>
840  void readOilWaterEffectiveParameters_(Container& dest,
841  const EclipseState& eclState,
842  unsigned satRegionIdx)
843  {
844  if (!hasOil || !hasWater)
845  // we don't read anything if either the water or the oil phase is not active
846  return;
847 
848  dest[satRegionIdx] = std::make_shared<OilWaterEffectiveTwoPhaseParams>();
849 
850  const auto tolcrit = eclState.runspec().saturationFunctionControls()
851  .minimumRelpermMobilityThreshold();
852 
853  const auto& tableManager = eclState.getTableManager();
854  auto& effParams = *dest[satRegionIdx];
855 
856  switch (eclState.runspec().saturationFunctionControls().family()) {
857  case SatFuncControls::KeywordFamily::Family_I:
858  {
859  const auto& swofTable = tableManager.getSwofTables().getTable<SwofTable>(satRegionIdx);
860  const std::vector<double> SwColumn = swofTable.getColumn("SW").vectorCopy();
861 
862  effParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn("KRW")));
863  effParams.setKrnSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn("KROW")));
864  effParams.setPcnwSamples(SwColumn, swofTable.getColumn("PCOW").vectorCopy());
865  effParams.finalize();
866  break;
867  }
868 
869  case SatFuncControls::KeywordFamily::Family_II:
870  {
871  const auto& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>(satRegionIdx);
872  const std::vector<double> SwColumn = swfnTable.getColumn("SW").vectorCopy();
873  effParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW")));
874  effParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy());
875 
876  if (!hasGas) {
877  const auto& sof2Table = tableManager.getSof2Tables().getTable<Sof2Table>(satRegionIdx);
878  // convert the saturations of the SOF2 keyword from oil to water saturations
879  std::vector<double> SwSamples(sof2Table.numRows());
880  for (size_t sampleIdx = 0; sampleIdx < sof2Table.numRows(); ++ sampleIdx)
881  SwSamples[sampleIdx] = 1 - sof2Table.get("SO", sampleIdx);
882 
883  effParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof2Table.getColumn("KRO")));
884  } else {
885  const auto& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>(satRegionIdx);
886  // convert the saturations of the SOF3 keyword from oil to water saturations
887  std::vector<double> SwSamples(sof3Table.numRows());
888  for (size_t sampleIdx = 0; sampleIdx < sof3Table.numRows(); ++ sampleIdx)
889  SwSamples[sampleIdx] = 1 - sof3Table.get("SO", sampleIdx);
890 
891  effParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof3Table.getColumn("KROW")));
892  }
893  effParams.finalize();
894  break;
895  }
896 
897  case SatFuncControls::KeywordFamily::Undefined:
898  throw std::domain_error("No valid saturation keyword family specified");
899  }
900  }
901 
902  template <class Container>
903  void readGasWaterEffectiveParameters_(Container& dest,
904  const EclipseState& eclState,
905  unsigned satRegionIdx)
906  {
907  if (!hasGas || !hasWater || hasOil)
908  // we don't read anything if either the gas or the water phase is not active or if oil is present
909  return;
910 
911  dest[satRegionIdx] = std::make_shared<GasWaterEffectiveTwoPhaseParams>();
912 
913  auto& effParams = *dest[satRegionIdx];
914 
915  const auto tolcrit = eclState.runspec().saturationFunctionControls()
916  .minimumRelpermMobilityThreshold();
917 
918  const auto& tableManager = eclState.getTableManager();
919 
920  switch (eclState.runspec().saturationFunctionControls().family()) {
921  case SatFuncControls::KeywordFamily::Family_I:
922  {
923  throw std::domain_error("Saturation keyword family I is not applicable for a gas-water system");
924  }
925 
926  case SatFuncControls::KeywordFamily::Family_II:
927  {
928  //Todo: allow also for Sgwfn table input as alternative to Sgfn and Swfn table input
929  const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable<SgfnTable>( satRegionIdx );
930  const SwfnTable& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>( satRegionIdx );
931 
932  std::vector<double> SwColumn = swfnTable.getColumn("SW").vectorCopy();
933 
934  effParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn("KRW")));
935  std::vector<double> SwSamples(sgfnTable.numRows());
936  for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx)
937  SwSamples[sampleIdx] = 1 - sgfnTable.get("SG", sampleIdx);
938  effParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn("KRG")));
939  //Capillary pressure is read from SWFN.
940  //For gas-water system the capillary pressure column values are set to 0 in SGFN
941  effParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy());
942  effParams.finalize();
943 
944  break;
945  }
946 
947  case SatFuncControls::KeywordFamily::Undefined:
948  throw std::domain_error("No valid saturation keyword family specified");
949  }
950  }
951 
952  template <class Container>
953  void readGasOilUnscaledPoints_(Container& dest,
954  std::shared_ptr<EclEpsConfig> config,
955  const EclipseState& /* eclState */,
956  unsigned satRegionIdx)
957  {
958  if (!hasGas || !hasOil)
959  // we don't read anything if either the gas or the oil phase is not active
960  return;
961 
962  dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
963  dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], *config, EclGasOilSystem);
964  }
965 
966  template <class Container>
967  void readOilWaterUnscaledPoints_(Container& dest,
968  std::shared_ptr<EclEpsConfig> config,
969  const EclipseState& /* eclState */,
970  unsigned satRegionIdx)
971  {
972  if (!hasOil || !hasWater)
973  // we don't read anything if either the water or the oil phase is not active
974  return;
975 
976  dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
977  dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], *config, EclOilWaterSystem);
978  }
979 
980  template <class Container>
981  void readGasWaterUnscaledPoints_(Container& dest,
982  std::shared_ptr<EclEpsConfig> config,
983  const EclipseState& /* eclState */,
984  unsigned satRegionIdx)
985  {
986  if (hasOil)
987  // we don't read anything if oil phase is active
988  return;
989 
990  dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
991  dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], *config, EclGasWaterSystem);
992  }
993 
994  std::tuple<EclEpsScalingPointsInfo<Scalar>,
995  EclEpsScalingPoints<Scalar>>
996  readScaledPoints_(const EclEpsConfig& config,
997  const EclipseState& eclState,
998  const EclEpsGridProperties& epsGridProperties,
999  unsigned elemIdx,
1000  EclTwoPhaseSystemType type)
1001  {
1002  unsigned satRegionIdx = epsGridProperties.satRegion( elemIdx );
1003 
1004  EclEpsScalingPointsInfo<Scalar> destInfo(unscaledEpsInfo_[satRegionIdx]);
1005  destInfo.extractScaled(eclState, epsGridProperties, elemIdx);
1006 
1007  EclEpsScalingPoints<Scalar> destPoint;
1008  destPoint.init(destInfo, config, type);
1009 
1010  return {destInfo, destPoint};
1011  }
1012 
1013  void initThreePhaseParams_(const EclipseState& /* eclState */,
1014  MaterialLawParams& materialParams,
1015  unsigned satRegionIdx,
1016  const EclEpsScalingPointsInfo<Scalar>& epsInfo,
1017  std::shared_ptr<OilWaterTwoPhaseHystParams> oilWaterParams,
1018  std::shared_ptr<GasOilTwoPhaseHystParams> gasOilParams,
1019  std::shared_ptr<GasWaterTwoPhaseHystParams> gasWaterParams)
1020  {
1021  materialParams.setApproach(threePhaseApproach_);
1022 
1023  switch (materialParams.approach()) {
1024  case EclMultiplexerApproach::EclStone1Approach: {
1025  auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclStone1Approach>();
1026  realParams.setGasOilParams(gasOilParams);
1027  realParams.setOilWaterParams(oilWaterParams);
1028  realParams.setSwl(epsInfo.Swl);
1029 
1030  if (!stoneEtas.empty()) {
1031  realParams.setEta(stoneEtas[satRegionIdx]);
1032  }
1033  else
1034  realParams.setEta(1.0);
1035  realParams.finalize();
1036  break;
1037  }
1038 
1039  case EclMultiplexerApproach::EclStone2Approach: {
1040  auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclStone2Approach>();
1041  realParams.setGasOilParams(gasOilParams);
1042  realParams.setOilWaterParams(oilWaterParams);
1043  realParams.setSwl(epsInfo.Swl);
1044  realParams.finalize();
1045  break;
1046  }
1047 
1048  case EclMultiplexerApproach::EclDefaultApproach: {
1049  auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclDefaultApproach>();
1050  realParams.setGasOilParams(gasOilParams);
1051  realParams.setOilWaterParams(oilWaterParams);
1052  realParams.setSwl(epsInfo.Swl);
1053  realParams.finalize();
1054  break;
1055  }
1056 
1057  case EclMultiplexerApproach::EclTwoPhaseApproach: {
1058  auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclTwoPhaseApproach>();
1059  realParams.setGasOilParams(gasOilParams);
1060  realParams.setOilWaterParams(oilWaterParams);
1061  realParams.setGasWaterParams(gasWaterParams);
1062  realParams.setApproach(twoPhaseApproach_);
1063  realParams.finalize();
1064  break;
1065  }
1066 
1067  case EclMultiplexerApproach::EclOnePhaseApproach: {
1068  // Nothing to do, no parameters.
1069  break;
1070  }
1071  }
1072  }
1073 
1074  // Relative permeability values not strictly greater than 'tolcrit' treated as zero.
1075  std::vector<double> normalizeKrValues_(const double tolcrit,
1076  const TableColumn& krValues) const
1077  {
1078  auto kr = krValues.vectorCopy();
1079  std::transform(kr.begin(), kr.end(), kr.begin(),
1080  [tolcrit](const double kri)
1081  {
1082  return (kri > tolcrit) ? kri : 0.0;
1083  });
1084 
1085  return kr;
1086  }
1087 
1088  bool enableEndPointScaling_;
1089  std::shared_ptr<EclHysteresisConfig> hysteresisConfig_;
1090 
1091  std::shared_ptr<EclEpsConfig> oilWaterEclEpsConfig_;
1092  std::vector<EclEpsScalingPointsInfo<Scalar>> unscaledEpsInfo_;
1093  OilWaterScalingInfoVector oilWaterScaledEpsInfoDrainage_;
1094 
1095  std::shared_ptr<EclEpsConfig> gasWaterEclEpsConfig_;
1096 
1097  GasOilScalingPointsVector gasOilUnscaledPointsVector_;
1098  OilWaterScalingPointsVector oilWaterUnscaledPointsVector_;
1099  GasWaterScalingPointsVector gasWaterUnscaledPointsVector_;
1100 
1101  GasOilEffectiveParamVector gasOilEffectiveParamVector_;
1102  OilWaterEffectiveParamVector oilWaterEffectiveParamVector_;
1103  GasWaterEffectiveParamVector gasWaterEffectiveParamVector_;
1104 
1105  EclMultiplexerApproach threePhaseApproach_ = EclMultiplexerApproach::EclDefaultApproach;
1106  // this attribute only makes sense for twophase simulations!
1107  enum EclTwoPhaseApproach twoPhaseApproach_ = EclTwoPhaseApproach::EclTwoPhaseGasOil;
1108 
1109  std::vector<MaterialLawParams> materialLawParams_;
1110 
1111  std::vector<int> satnumRegionArray_;
1112  std::vector<int> imbnumRegionArray_;
1113  std::vector<Scalar> stoneEtas;
1114 
1115  bool hasGas;
1116  bool hasOil;
1117  bool hasWater;
1118 
1119  std::shared_ptr<EclEpsConfig> gasOilConfig;
1120  std::shared_ptr<EclEpsConfig> oilWaterConfig;
1121  std::shared_ptr<EclEpsConfig> gasWaterConfig;
1122 };
1123 } // namespace Opm
1124 
1125 #endif
Specifies the configuration used by the endpoint scaling code.
This material law takes a material law defined for unscaled saturation and converts it to a material ...
Specifies the configuration used by the ECL kr/pC hysteresis code.
This material law implements the hysteresis model of the ECL file format.
Implements a multiplexer class that provides all three phase capillary pressure laws used by the ECLi...
Implementation for the parameters required by the material law for two-phase simulations.
This file contains helper classes for the material laws.
Implementation of a tabulated, piecewise linear capillary pressure law.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: EclEpsGridProperties.hpp:76
This material law takes a material law defined for unscaled saturation and converts it to a material ...
Definition: EclEpsTwoPhaseLaw.hpp:53
This material law implements the hysteresis model of the ECL file format.
Definition: EclHysteresisTwoPhaseLaw.hpp:41
Provides an simple way to create and manage the material law objects for a complete ECL deck.
Definition: EclMaterialLawManager.hpp:69
const MaterialLawParams & connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const
Returns a material parameter object for a given element and saturation region.
Definition: EclMaterialLawManager.hpp:484
Scalar applySwatinit(unsigned elemIdx, Scalar pcow, Scalar Sw)
Modify the initial condition according to the SWATINIT keyword.
Definition: EclMaterialLawManager.hpp:410
Implements a multiplexer class that provides all three phase capillary pressure laws used by the ECLi...
Definition: EclMultiplexerMaterial.hpp:58
static void capillaryPressures(ContainerT &values, const Params &params, const FluidState &fluidState)
Implements the multiplexer three phase capillary pressure law used by the ECLipse simulator.
Definition: EclMultiplexerMaterial.hpp:135
static void updateHysteresis(Params &params, const FluidState &fluidState)
Update the hysteresis parameters after a time step.
Definition: EclMultiplexerMaterial.hpp:545
Implementation of a tabulated, piecewise linear capillary pressure law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:51
ParamsT Params
The type of the parameter objects for this law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:59
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: SimpleModularFluidState.hpp:104
A generic traits class for two-phase material laws.
Definition: MaterialTraits.hpp:60