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