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