28 #error "The opm-parser module is required to use the ECL material manager!" 31 #ifndef OPM_ECL_MATERIAL_LAW_MANAGER_HPP 32 #define OPM_ECL_MATERIAL_LAW_MANAGER_HPP 45 #include <opm/common/Exceptions.hpp> 46 #include <opm/common/ErrorMacros.hpp> 48 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp> 49 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp> 50 #include <opm/parser/eclipse/Deck/Deck.hpp> 63 template <
class TraitsT>
67 typedef TraitsT Traits;
68 typedef typename Traits::Scalar Scalar;
69 enum { waterPhaseIdx = Traits::wettingPhaseIdx };
70 enum { oilPhaseIdx = Traits::nonWettingPhaseIdx };
71 enum { gasPhaseIdx = Traits::gasPhaseIdx };
72 enum { numPhases = Traits::numPhases };
86 typedef typename GasOilEpsTwoPhaseLaw::Params GasOilEpsTwoPhaseParams;
87 typedef typename OilWaterEpsTwoPhaseLaw::Params OilWaterEpsTwoPhaseParams;
92 typedef typename GasOilTwoPhaseLaw::Params GasOilTwoPhaseHystParams;
93 typedef typename OilWaterTwoPhaseLaw::Params OilWaterTwoPhaseHystParams;
98 typedef typename MaterialLaw::Params MaterialLawParams;
102 typedef std::vector<std::shared_ptr<GasOilEffectiveTwoPhaseParams> > GasOilEffectiveParamVector;
103 typedef std::vector<std::shared_ptr<OilWaterEffectiveTwoPhaseParams> > OilWaterEffectiveParamVector;
104 typedef std::vector<std::shared_ptr<EclEpsScalingPoints<Scalar> > > GasOilScalingPointsVector;
105 typedef std::vector<std::shared_ptr<EclEpsScalingPoints<Scalar> > > OilWaterScalingPointsVector;
106 typedef std::vector<std::shared_ptr<EclEpsScalingPointsInfo<Scalar> > > GasOilScalingInfoVector;
107 typedef std::vector<std::shared_ptr<EclEpsScalingPointsInfo<Scalar> > > OilWaterScalingInfoVector;
108 typedef std::vector<std::shared_ptr<GasOilTwoPhaseHystParams> > GasOilParamVector;
109 typedef std::vector<std::shared_ptr<OilWaterTwoPhaseHystParams> > OilWaterParamVector;
110 typedef std::vector<std::shared_ptr<MaterialLawParams> > MaterialLawParamsVector;
116 void initFromDeck(
const Opm::Deck& deck,
117 const Opm::EclipseState& eclState,
118 const std::vector<int>& compressedToCartesianElemIdx)
121 const size_t numSatRegions = eclState.runspec().tabdims().getNumSatTables();
122 size_t numCompressedElems = compressedToCartesianElemIdx.size();
126 satnumRegionArray_.resize(numCompressedElems);
127 if (eclState.get3DProperties().hasDeckIntGridProperty(
"SATNUM")) {
128 const auto& satnumRawData = eclState.get3DProperties().getIntGridProperty(
"SATNUM").getData();
129 for (
unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
130 unsigned cartesianElemIdx =
static_cast<unsigned>(compressedToCartesianElemIdx[elemIdx]);
131 satnumRegionArray_[elemIdx] = satnumRawData[cartesianElemIdx] - 1;
135 std::fill(satnumRegionArray_.begin(), satnumRegionArray_.end(), 0);
137 readGlobalEpsOptions_(deck, eclState);
138 readGlobalHysteresisOptions_(deck);
139 readGlobalThreePhaseOptions_(deck);
141 unscaledEpsInfo_.resize(numSatRegions);
142 for (
unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx)
143 unscaledEpsInfo_[satRegionIdx].extractUnscaled(deck, eclState, satRegionIdx);
145 initParamsForElements_(deck, eclState, compressedToCartesianElemIdx, satnumRegionArray_);
160 auto& elemScaledEpsInfo = *oilWaterScaledEpsInfoDrainage_[elemIdx];
165 Sw = elemScaledEpsInfo.Swu;
168 if (Sw <= elemScaledEpsInfo.Swl)
169 Sw = elemScaledEpsInfo.Swl;
185 fs.setSaturation(waterPhaseIdx, Sw);
186 fs.setSaturation(gasPhaseIdx, 0);
187 fs.setSaturation(oilPhaseIdx, 0);
188 Scalar pc[numPhases] = { 0 };
191 Scalar pcowAtSw = pc[oilPhaseIdx] - pc[waterPhaseIdx];
192 if (pcowAtSw > 0.0) {
193 elemScaledEpsInfo.maxPcow *= pcow/pcowAtSw;
194 auto& elemEclEpsScalingPoints = oilWaterScaledEpsPointsDrainage(elemIdx);
195 elemEclEpsScalingPoints.init(elemScaledEpsInfo, *oilWaterEclEpsConfig_, Opm::EclOilWaterSystem);
202 bool enableEndPointScaling()
const 203 {
return enableEndPointScaling_; }
205 bool enableHysteresis()
const 206 {
return hysteresisConfig_->enableHysteresis(); }
208 MaterialLawParams& materialLawParams(
unsigned elemIdx)
210 assert(0 <= elemIdx && elemIdx < materialLawParams_.size());
211 return *materialLawParams_[elemIdx];
214 const MaterialLawParams& materialLawParams(
unsigned elemIdx)
const 216 assert(0 <= elemIdx && elemIdx < materialLawParams_.size());
217 return *materialLawParams_[elemIdx];
230 MaterialLawParams& mlp = *materialLawParams_[elemIdx];
232 if (enableHysteresis()) {
233 OPM_MESSAGE(
"Warning: Using non-defautl satnum regions for conenction is not tested in combination with hysteresis");
239 switch (mlp.approach()) {
240 case EclStone1Approach: {
241 auto& realParams = mlp.template getRealParams<Opm::EclStone1Approach>();
243 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
244 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
245 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
246 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
256 case EclStone2Approach: {
257 auto& realParams = mlp.template getRealParams<Opm::EclStone2Approach>();
258 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
259 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
260 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
261 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
271 case EclDefaultApproach: {
272 auto& realParams = mlp.template getRealParams<Opm::EclDefaultApproach>();
273 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
274 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
275 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
276 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
286 case EclTwoPhaseApproach: {
287 auto& realParams = mlp.template getRealParams<Opm::EclTwoPhaseApproach>();
288 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
289 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
290 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
291 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
302 OPM_THROW(std::logic_error,
"Enum value for material approach unknown!");
308 int satnumRegionIdx(
unsigned elemIdx)
const {
309 return satnumRegionArray_[elemIdx];
312 std::shared_ptr<MaterialLawParams>& materialLawParamsPointerReferenceHack(
unsigned elemIdx)
314 assert(0 <= elemIdx && elemIdx < materialLawParams_.size());
315 return materialLawParams_[elemIdx];
318 template <
class Flu
idState>
319 void updateHysteresis(
const FluidState& fluidState,
unsigned elemIdx)
321 if (!enableHysteresis())
324 auto threePhaseParams = materialLawParams_[elemIdx];
328 void oilWaterHysteresisParams(Scalar& pcSwMdc,
330 unsigned elemIdx)
const 332 if (!enableHysteresis()) {
333 OPM_THROW(std::runtime_error,
"Cannot get hysteresis parameters if hysteresis not enabled.");
335 const auto& params = materialLawParams(elemIdx);
336 MaterialLaw::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, params);
339 void setOilWaterHysteresisParams(
const Scalar& pcSwMdc,
340 const Scalar& krnSwMdc,
343 if (!enableHysteresis()) {
344 OPM_THROW(std::runtime_error,
"Cannot set hysteresis parameters if hysteresis not enabled.");
346 auto& params = materialLawParams(elemIdx);
347 MaterialLaw::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, params);
350 void gasOilHysteresisParams(Scalar& pcSwMdc,
352 unsigned elemIdx)
const 354 if (!enableHysteresis()) {
355 OPM_THROW(std::runtime_error,
"Cannot get hysteresis parameters if hysteresis not enabled.");
357 const auto& params = materialLawParams(elemIdx);
358 MaterialLaw::gasOilHysteresisParams(pcSwMdc, krnSwMdc, params);
361 void setGasOilHysteresisParams(
const Scalar& pcSwMdc,
362 const Scalar& krnSwMdc,
365 if (!enableHysteresis()) {
366 OPM_THROW(std::runtime_error,
"Cannot set hysteresis parameters if hysteresis not enabled.");
368 auto& params = materialLawParams(elemIdx);
369 MaterialLaw::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, params);
372 EclEpsScalingPoints<Scalar>& oilWaterScaledEpsPointsDrainage(
unsigned elemIdx)
374 auto& materialParams = *materialLawParams_[elemIdx];
375 switch (materialParams.approach()) {
376 case EclStone1Approach: {
377 auto& realParams = materialParams.template getRealParams<Opm::EclStone1Approach>();
378 return realParams.oilWaterParams().drainageParams().scaledPoints();
381 case EclStone2Approach: {
382 auto& realParams = materialParams.template getRealParams<Opm::EclStone2Approach>();
383 return realParams.oilWaterParams().drainageParams().scaledPoints();
386 case EclDefaultApproach: {
387 auto& realParams = materialParams.template getRealParams<Opm::EclDefaultApproach>();
388 return realParams.oilWaterParams().drainageParams().scaledPoints();
391 case EclTwoPhaseApproach: {
392 auto& realParams = materialParams.template getRealParams<Opm::EclTwoPhaseApproach>();
393 return realParams.oilWaterParams().drainageParams().scaledPoints();
396 OPM_THROW(std::logic_error,
"Enum value for material approach unknown!");
402 return *oilWaterScaledEpsInfoDrainage_[elemIdx];
405 std::shared_ptr<EclEpsScalingPointsInfo<Scalar> >& oilWaterScaledEpsInfoDrainagePointerReferenceHack(
unsigned elemIdx)
407 return oilWaterScaledEpsInfoDrainage_[elemIdx];
410 void readGlobalEpsOptions_(
const Opm::Deck& deck,
const Opm::EclipseState& eclState)
412 oilWaterEclEpsConfig_ = std::make_shared<Opm::EclEpsConfig>();
413 oilWaterEclEpsConfig_-> initFromDeck(deck, eclState, Opm::EclOilWaterSystem);
415 enableEndPointScaling_ = deck.hasKeyword(
"ENDSCALE");
418 void readGlobalHysteresisOptions_(
const Opm::Deck& deck)
420 hysteresisConfig_ = std::make_shared<Opm::EclHysteresisConfig>();
421 hysteresisConfig_->initFromDeck(deck);
424 void readGlobalThreePhaseOptions_(
const Opm::Deck& deck)
426 bool gasEnabled = deck.hasKeyword(
"GAS");
427 bool oilEnabled = deck.hasKeyword(
"OIL");
428 bool waterEnabled = deck.hasKeyword(
"WATER");
433 + (waterEnabled?1:0);
436 OPM_THROW(std::runtime_error,
437 "At least two fluid phases must be enabled. (Is: " << numEnabled <<
")");
439 if (numEnabled == 2) {
440 threePhaseApproach_ = Opm::EclTwoPhaseApproach;
442 twoPhaseApproach_ = Opm::EclTwoPhaseOilWater;
443 else if (!oilEnabled)
444 twoPhaseApproach_ = Opm::EclTwoPhaseGasWater;
445 else if (!waterEnabled)
446 twoPhaseApproach_ = Opm::EclTwoPhaseGasOil;
449 assert(numEnabled == 3);
451 threePhaseApproach_ = Opm::EclDefaultApproach;
452 if (deck.hasKeyword(
"STONE") || deck.hasKeyword(
"STONE2"))
453 threePhaseApproach_ = Opm::EclStone2Approach;
454 else if (deck.hasKeyword(
"STONE1"))
455 threePhaseApproach_ = Opm::EclStone1Approach;
459 void initParamsForElements_(
const Deck& deck,
const EclipseState& eclState,
460 const std::vector<int>& compressedToCartesianElemIdx,
461 const std::vector<int>& satnumRegionArray)
463 const size_t numSatRegions = eclState.runspec().tabdims().getNumSatTables();
464 unsigned numCompressedElems =
static_cast<unsigned>(compressedToCartesianElemIdx.size());
468 auto gasOilConfig = std::make_shared<Opm::EclEpsConfig>();
469 auto oilWaterConfig = std::make_shared<Opm::EclEpsConfig>();
470 gasOilConfig->initFromDeck(deck, eclState, Opm::EclGasOilSystem);
471 oilWaterConfig->initFromDeck(deck, eclState, Opm::EclOilWaterSystem);
474 gasOilUnscaledPointsVector_.resize(numSatRegions);
475 oilWaterUnscaledPointsVector_.resize(numSatRegions);
476 gasOilEffectiveParamVector_.resize(numSatRegions);
477 oilWaterEffectiveParamVector_.resize(numSatRegions);
478 for (
unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) {
480 readGasOilUnscaledPoints_(gasOilUnscaledPointsVector_, gasOilConfig, deck, eclState, satRegionIdx);
481 readOilWaterUnscaledPoints_(oilWaterUnscaledPointsVector_, oilWaterConfig, deck, eclState, satRegionIdx);
484 readGasOilEffectiveParameters_(gasOilEffectiveParamVector_, deck, eclState, satRegionIdx);
485 readOilWaterEffectiveParameters_(oilWaterEffectiveParamVector_, deck, eclState, satRegionIdx);
488 unscaledEpsInfo_[satRegionIdx].extractUnscaled(deck, eclState, satRegionIdx);
494 GasOilScalingInfoVector gasOilScaledInfoVector(numCompressedElems);
495 oilWaterScaledEpsInfoDrainage_.resize(numCompressedElems);
496 GasOilScalingInfoVector gasOilScaledImbInfoVector;
497 OilWaterScalingInfoVector oilWaterScaledImbInfoVector;
499 GasOilScalingPointsVector gasOilScaledPointsVector(numCompressedElems);
500 GasOilScalingPointsVector oilWaterScaledEpsPointsDrainage(numCompressedElems);
501 GasOilScalingPointsVector gasOilScaledImbPointsVector;
502 OilWaterScalingPointsVector oilWaterScaledImbPointsVector;
504 if (enableHysteresis()) {
505 gasOilScaledImbInfoVector.resize(numCompressedElems);
506 gasOilScaledImbPointsVector.resize(numCompressedElems);
507 oilWaterScaledImbInfoVector.resize(numCompressedElems);
508 oilWaterScaledImbPointsVector.resize(numCompressedElems);
511 EclEpsGridProperties epsGridProperties, epsImbGridProperties;
512 epsGridProperties.initFromDeck(deck, eclState,
false);
513 if (enableHysteresis())
514 epsImbGridProperties.initFromDeck(deck, eclState,
true);
515 for (
unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
516 unsigned cartElemIdx =
static_cast<unsigned>(compressedToCartesianElemIdx[elemIdx]);
517 readGasOilScaledPoints_(gasOilScaledInfoVector,
518 gasOilScaledPointsVector,
524 readOilWaterScaledPoints_(oilWaterScaledEpsInfoDrainage_,
525 oilWaterScaledEpsPointsDrainage,
532 if (enableHysteresis()) {
533 readGasOilScaledPoints_(gasOilScaledImbInfoVector,
534 gasOilScaledImbPointsVector,
537 epsImbGridProperties,
540 readOilWaterScaledPoints_(oilWaterScaledImbInfoVector,
541 oilWaterScaledImbPointsVector,
544 epsImbGridProperties,
551 GasOilParamVector gasOilParams(numCompressedElems);
552 OilWaterParamVector oilWaterParams(numCompressedElems);
553 GasOilParamVector gasOilImbParams;
554 OilWaterParamVector oilWaterImbParams;
556 if (enableHysteresis()) {
557 gasOilImbParams.resize(numCompressedElems);
558 oilWaterImbParams.resize(numCompressedElems);
561 bool hasGas = deck.hasKeyword(
"GAS");
562 bool hasOil = deck.hasKeyword(
"OIL");
563 bool hasWater = deck.hasKeyword(
"WATER");
565 const auto& imbnumData = eclState.get3DProperties().getIntGridProperty(
"IMBNUM").getData();
566 assert(numCompressedElems == satnumRegionArray.size());
567 for (
unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
568 unsigned satRegionIdx =
static_cast<unsigned>(satnumRegionArray[elemIdx]);
570 gasOilParams[elemIdx] = std::make_shared<GasOilTwoPhaseHystParams>();
571 oilWaterParams[elemIdx] = std::make_shared<OilWaterTwoPhaseHystParams>();
573 gasOilParams[elemIdx]->setConfig(hysteresisConfig_);
574 oilWaterParams[elemIdx]->setConfig(hysteresisConfig_);
576 if (hasGas && hasOil) {
577 auto gasOilDrainParams = std::make_shared<GasOilEpsTwoPhaseParams>();
578 gasOilDrainParams->setConfig(gasOilConfig);
579 gasOilDrainParams->setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
580 gasOilDrainParams->setScaledPoints(gasOilScaledPointsVector[elemIdx]);
581 gasOilDrainParams->setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
582 gasOilDrainParams->finalize();
584 gasOilParams[elemIdx]->setDrainageParams(gasOilDrainParams,
585 *gasOilScaledInfoVector[elemIdx],
589 if (hasOil && hasWater) {
590 auto oilWaterDrainParams = std::make_shared<OilWaterEpsTwoPhaseParams>();
591 oilWaterDrainParams->setConfig(oilWaterConfig);
592 oilWaterDrainParams->setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
593 oilWaterDrainParams->setScaledPoints(oilWaterScaledEpsPointsDrainage[elemIdx]);
594 oilWaterDrainParams->setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
595 oilWaterDrainParams->finalize();
597 oilWaterParams[elemIdx]->setDrainageParams(oilWaterDrainParams,
598 *oilWaterScaledEpsInfoDrainage_[elemIdx],
602 if (enableHysteresis()) {
603 unsigned imbRegionIdx =
static_cast<unsigned>(imbnumData[elemIdx]) - 1;
605 if (hasGas && hasOil) {
606 auto gasOilImbParamsHyst = std::make_shared<GasOilEpsTwoPhaseParams>();
607 gasOilImbParamsHyst->setConfig(gasOilConfig);
608 gasOilImbParamsHyst->setUnscaledPoints(gasOilUnscaledPointsVector_[imbRegionIdx]);
609 gasOilImbParamsHyst->setScaledPoints(gasOilScaledImbPointsVector[elemIdx]);
610 gasOilImbParamsHyst->setEffectiveLawParams(gasOilEffectiveParamVector_[imbRegionIdx]);
611 gasOilImbParamsHyst->finalize();
613 gasOilParams[elemIdx]->setImbibitionParams(gasOilImbParamsHyst,
614 *gasOilScaledImbInfoVector[elemIdx],
618 if (hasOil && hasWater) {
619 auto oilWaterImbParamsHyst = std::make_shared<OilWaterEpsTwoPhaseParams>();
620 oilWaterImbParamsHyst->setConfig(oilWaterConfig);
621 oilWaterImbParamsHyst->setUnscaledPoints(oilWaterUnscaledPointsVector_[imbRegionIdx]);
622 oilWaterImbParamsHyst->setScaledPoints(oilWaterScaledImbPointsVector[elemIdx]);
623 oilWaterImbParamsHyst->setEffectiveLawParams(oilWaterEffectiveParamVector_[imbRegionIdx]);
624 oilWaterImbParamsHyst->finalize();
626 oilWaterParams[elemIdx]->setImbibitionParams(oilWaterImbParamsHyst,
627 *gasOilScaledImbInfoVector[elemIdx],
632 if (hasGas && hasOil)
633 gasOilParams[elemIdx]->finalize();
635 if (hasOil && hasWater)
636 oilWaterParams[elemIdx]->finalize();
640 materialLawParams_.resize(numCompressedElems);
641 for (
unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
642 materialLawParams_[elemIdx] = std::make_shared<MaterialLawParams>();
643 unsigned satRegionIdx =
static_cast<unsigned>(satnumRegionArray[elemIdx]);
645 initThreePhaseParams_(deck,
647 *materialLawParams_[elemIdx],
649 *oilWaterScaledEpsInfoDrainage_[elemIdx],
650 oilWaterParams[elemIdx],
651 gasOilParams[elemIdx]);
653 materialLawParams_[elemIdx]->finalize();
661 enum SaturationFunctionFamily {
667 SaturationFunctionFamily getSaturationFunctionFamily(
const Opm::Deck& deck,
const Opm::EclipseState& eclState)
const 669 const auto& tableManager = eclState.getTableManager();
670 const TableContainer& swofTables = tableManager.getSwofTables();
671 const TableContainer& slgofTables= tableManager.getSlgofTables();
672 const TableContainer& sgofTables = tableManager.getSgofTables();
673 const TableContainer& swfnTables = tableManager.getSwfnTables();
674 const TableContainer& sgfnTables = tableManager.getSgfnTables();
675 const TableContainer& sof3Tables = tableManager.getSof3Tables();
677 bool hasGas = deck.hasKeyword(
"GAS");
678 bool hasOil = deck.hasKeyword(
"OIL");
679 bool hasWater = deck.hasKeyword(
"WATER");
681 bool family1 =
false;
682 bool family2 =
false;
685 family1 = !swofTables.empty();
687 throw std::runtime_error(
"only SWOF is supporeted for oil-water two-phase simulations");
690 else if (!hasWater) {
692 family1 = !sgofTables.empty();
694 throw std::runtime_error(
"only SGOF is supporeted for oil-gas two-phase simulations");
699 throw std::runtime_error(
"water-gas two-phase simulations are currently not supported");
703 family1 = (!sgofTables.empty() || !slgofTables.empty()) && !swofTables.empty();
704 family2 = !swfnTables.empty() && !sgfnTables.empty() && !sof3Tables.empty();
707 if (family1 && family2) {
708 throw std::invalid_argument(
"Saturation families should not be mixed \n" 709 "Use either SGOF and SWOF or SGFN, SWFN and SOF3");
712 if (!family1 && !family2) {
713 throw std::invalid_argument(
"Saturations function must be specified using either " 714 "family 1 or family 2 keywords \n" 715 "Use either SGOF and SWOF or SGFN, SWFN and SOF3" );
718 if (family1 && !family2)
719 return SaturationFunctionFamily::FamilyI;
720 else if (family2 && !family1)
721 return SaturationFunctionFamily::FamilyII;
722 return SaturationFunctionFamily::noFamily;
725 template <
class Container>
726 void readGasOilEffectiveParameters_(Container& dest,
727 const Opm::Deck& deck,
728 const Opm::EclipseState& eclState,
729 unsigned satRegionIdx)
731 bool hasGas = deck.hasKeyword(
"GAS");
732 bool hasOil = deck.hasKeyword(
"OIL");
734 if (!hasGas || !hasOil)
738 dest[satRegionIdx] = std::make_shared<GasOilEffectiveTwoPhaseParams>();
740 auto& effParams = *dest[satRegionIdx];
744 Scalar Swco = unscaledEpsInfo_[satRegionIdx].Swl;
745 const auto& tableManager = eclState.getTableManager();
747 switch (getSaturationFunctionFamily(deck, eclState)) {
750 const TableContainer& sgofTables = tableManager.getSgofTables();
751 const TableContainer& slgofTables = tableManager.getSlgofTables();
752 if (!sgofTables.empty())
753 readGasOilEffectiveParametersSgof_(effParams,
755 sgofTables.getTable<SgofTable>(satRegionIdx));
756 else if (!slgofTables.empty())
757 readGasOilEffectiveParametersSlgof_(effParams,
759 slgofTables.getTable<SlgofTable>(satRegionIdx));
765 const Sof3Table& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>( satRegionIdx );
766 const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable<SgfnTable>( satRegionIdx );
767 readGasOilEffectiveParametersFamily2_(effParams,
776 throw std::domain_error(
"No valid saturation keyword family specified");
781 void readGasOilEffectiveParametersSgof_(GasOilEffectiveTwoPhaseParams& effParams,
783 const Opm::SgofTable& sgofTable)
786 std::vector<double> SoSamples(sgofTable.numRows());
787 std::vector<double> SoKroSamples(sgofTable.numRows());
788 for (
size_t sampleIdx = 0; sampleIdx < sgofTable.numRows(); ++ sampleIdx) {
789 SoSamples[sampleIdx] = 1 - sgofTable.get(
"SG", sampleIdx);
790 SoKroSamples[sampleIdx] = SoSamples[sampleIdx] - Swco;
793 effParams.setKrwSamples(SoKroSamples, sgofTable.getColumn(
"KROG").vectorCopy());
794 effParams.setKrnSamples(SoSamples, sgofTable.getColumn(
"KRG").vectorCopy());
795 effParams.setPcnwSamples(SoSamples, sgofTable.getColumn(
"PCOG").vectorCopy());
796 effParams.finalize();
799 void readGasOilEffectiveParametersSlgof_(GasOilEffectiveTwoPhaseParams& effParams,
801 const Opm::SlgofTable& slgofTable)
804 std::vector<double> SoSamples(slgofTable.numRows());
805 std::vector<double> SoKroSamples(slgofTable.numRows());
806 for (
size_t sampleIdx = 0; sampleIdx < slgofTable.numRows(); ++ sampleIdx) {
807 SoSamples[sampleIdx] = slgofTable.get(
"SL", sampleIdx);
808 SoKroSamples[sampleIdx] = slgofTable.get(
"SL", sampleIdx) - Swco;
811 effParams.setKrwSamples(SoKroSamples, slgofTable.getColumn(
"KROG").vectorCopy());
812 effParams.setKrnSamples(SoSamples, slgofTable.getColumn(
"KRG").vectorCopy());
813 effParams.setPcnwSamples(SoSamples, slgofTable.getColumn(
"PCOG").vectorCopy());
814 effParams.finalize();
817 void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
819 const Opm::Sof3Table& sof3Table,
820 const Opm::SgfnTable& sgfnTable)
823 std::vector<double> SoSamples(sgfnTable.numRows());
824 std::vector<double> SoColumn = sof3Table.getColumn(
"SO").vectorCopy();
825 for (
size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) {
826 SoSamples[sampleIdx] = 1 - sgfnTable.get(
"SG", sampleIdx);
829 effParams.setKrwSamples(SoColumn, sof3Table.getColumn(
"KROG").vectorCopy());
830 effParams.setKrnSamples(SoSamples, sgfnTable.getColumn(
"KRG").vectorCopy());
831 effParams.setPcnwSamples(SoSamples, sgfnTable.getColumn(
"PCOG").vectorCopy());
832 effParams.finalize();
835 template <
class Container>
836 void readOilWaterEffectiveParameters_(Container& dest,
837 const Opm::Deck& deck,
838 const Opm::EclipseState& eclState,
839 unsigned satRegionIdx)
841 bool hasWater = deck.hasKeyword(
"WATER");
842 bool hasOil = deck.hasKeyword(
"OIL");
844 if (!hasOil || !hasWater)
848 dest[satRegionIdx] = std::make_shared<OilWaterEffectiveTwoPhaseParams>();
850 const auto& tableManager = eclState.getTableManager();
851 auto& effParams = *dest[satRegionIdx];
853 switch (getSaturationFunctionFamily(deck, eclState)) {
855 const auto& swofTable = tableManager.getSwofTables().getTable<SwofTable>(satRegionIdx);
856 std::vector<double> SwColumn = swofTable.getColumn(
"SW").vectorCopy();
858 effParams.setKrwSamples(SwColumn, swofTable.getColumn(
"KRW").vectorCopy());
859 effParams.setKrnSamples(SwColumn, swofTable.getColumn(
"KROW").vectorCopy());
860 effParams.setPcnwSamples(SwColumn, swofTable.getColumn(
"PCOW").vectorCopy());
861 effParams.finalize();
866 const auto& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>(satRegionIdx);
867 const auto& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>(satRegionIdx);
868 std::vector<double> SwColumn = swfnTable.getColumn(
"SW").vectorCopy();
871 std::vector<double> SwSamples(sof3Table.numRows());
872 for (
size_t sampleIdx = 0; sampleIdx < sof3Table.numRows(); ++ sampleIdx)
873 SwSamples[sampleIdx] = 1 - sof3Table.get(
"SO", sampleIdx);
875 effParams.setKrwSamples(SwColumn, swfnTable.getColumn(
"KRW").vectorCopy());
876 effParams.setKrnSamples(SwSamples, sof3Table.getColumn(
"KROW").vectorCopy());
877 effParams.setPcnwSamples(SwColumn, swfnTable.getColumn(
"PCOW").vectorCopy());
878 effParams.finalize();
884 throw std::domain_error(
"No valid saturation keyword family specified");
890 template <
class Container>
891 void readGasOilUnscaledPoints_(Container& dest,
892 std::shared_ptr<EclEpsConfig> config,
893 const Opm::Deck& deck,
894 const Opm::EclipseState& ,
895 unsigned satRegionIdx)
897 bool hasGas = deck.hasKeyword(
"GAS");
898 bool hasOil = deck.hasKeyword(
"OIL");
900 if (!hasGas || !hasOil)
904 dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
905 dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], *config, EclGasOilSystem);
908 template <
class Container>
909 void readOilWaterUnscaledPoints_(Container& dest,
910 std::shared_ptr<EclEpsConfig> config,
911 const Opm::Deck& deck,
912 const Opm::EclipseState& ,
913 unsigned satRegionIdx)
915 bool hasWater = deck.hasKeyword(
"WATER");
916 bool hasOil = deck.hasKeyword(
"OIL");
918 if (!hasOil || !hasWater)
922 dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
923 dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], *config, EclOilWaterSystem);
926 template <
class InfoContainer,
class Po
intsContainer>
927 void readGasOilScaledPoints_(InfoContainer& destInfo,
928 PointsContainer& destPoints,
929 std::shared_ptr<EclEpsConfig> config,
930 const Opm::EclipseState& eclState,
931 const EclEpsGridProperties& epsGridProperties,
933 unsigned cartElemIdx)
935 unsigned satRegionIdx =
static_cast<unsigned>((*epsGridProperties.satnum)[cartElemIdx]) - 1;
937 destInfo[elemIdx] = std::make_shared<EclEpsScalingPointsInfo<Scalar> >(unscaledEpsInfo_[satRegionIdx]);
938 destInfo[elemIdx]->extractScaled(eclState, epsGridProperties, cartElemIdx);
940 destPoints[elemIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
941 destPoints[elemIdx]->init(*destInfo[elemIdx], *config, EclGasOilSystem);
944 template <
class InfoContainer,
class Po
intsContainer>
945 void readOilWaterScaledPoints_(InfoContainer& destInfo,
946 PointsContainer& destPoints,
947 std::shared_ptr<EclEpsConfig> config,
948 const Opm::EclipseState& eclState,
949 const EclEpsGridProperties& epsGridProperties,
951 unsigned cartElemIdx)
953 unsigned satRegionIdx =
static_cast<unsigned>((*epsGridProperties.satnum)[cartElemIdx]) - 1;
955 destInfo[elemIdx] = std::make_shared<EclEpsScalingPointsInfo<Scalar> >(unscaledEpsInfo_[satRegionIdx]);
956 destInfo[elemIdx]->extractScaled(eclState, epsGridProperties, cartElemIdx);
958 destPoints[elemIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
959 destPoints[elemIdx]->init(*destInfo[elemIdx], *config, EclOilWaterSystem);
962 void initThreePhaseParams_(
const Opm::Deck& deck,
963 const Opm::EclipseState& ,
964 MaterialLawParams& materialParams,
965 unsigned satRegionIdx,
966 const EclEpsScalingPointsInfo<Scalar>& epsInfo,
967 std::shared_ptr<OilWaterTwoPhaseHystParams> oilWaterParams,
968 std::shared_ptr<GasOilTwoPhaseHystParams> gasOilParams)
970 materialParams.setApproach(threePhaseApproach_);
972 switch (materialParams.approach()) {
973 case EclStone1Approach: {
974 auto& realParams = materialParams.template getRealParams<Opm::EclStone1Approach>();
975 realParams.setGasOilParams(gasOilParams);
976 realParams.setOilWaterParams(oilWaterParams);
977 realParams.setSwl(epsInfo.Swl);
979 if (deck.hasKeyword(
"STONE1EX")) {
981 deck.getKeyword(
"STONE1EX").getRecord(satRegionIdx).getItem(0).getSIDouble(0);
982 realParams.setEta(eta);
985 realParams.setEta(1.0);
986 realParams.finalize();
990 case EclStone2Approach: {
991 auto& realParams = materialParams.template getRealParams<Opm::EclStone2Approach>();
992 realParams.setGasOilParams(gasOilParams);
993 realParams.setOilWaterParams(oilWaterParams);
994 realParams.setSwl(epsInfo.Swl);
995 realParams.finalize();
999 case EclDefaultApproach: {
1000 auto& realParams = materialParams.template getRealParams<Opm::EclDefaultApproach>();
1001 realParams.setGasOilParams(gasOilParams);
1002 realParams.setOilWaterParams(oilWaterParams);
1003 realParams.setSwl(epsInfo.Swl);
1004 realParams.finalize();
1008 case EclTwoPhaseApproach: {
1009 auto& realParams = materialParams.template getRealParams<Opm::EclTwoPhaseApproach>();
1010 realParams.setGasOilParams(gasOilParams);
1011 realParams.setOilWaterParams(oilWaterParams);
1012 realParams.setApproach(twoPhaseApproach_);
1013 realParams.finalize();
1019 bool enableEndPointScaling_;
1020 std::shared_ptr<EclHysteresisConfig> hysteresisConfig_;
1022 std::shared_ptr<EclEpsConfig> oilWaterEclEpsConfig_;
1023 std::vector<Opm::EclEpsScalingPointsInfo<Scalar>> unscaledEpsInfo_;
1024 OilWaterScalingInfoVector oilWaterScaledEpsInfoDrainage_;
1026 GasOilScalingPointsVector gasOilUnscaledPointsVector_;
1027 OilWaterScalingPointsVector oilWaterUnscaledPointsVector_;
1028 GasOilEffectiveParamVector gasOilEffectiveParamVector_;
1029 OilWaterEffectiveParamVector oilWaterEffectiveParamVector_;
1031 Opm::EclMultiplexerApproach threePhaseApproach_;
1033 enum EclTwoPhaseApproach twoPhaseApproach_;
1035 std::vector<std::shared_ptr<MaterialLawParams> > materialLawParams_;
1037 std::vector<int> satnumRegionArray_;
This material law implements the hysteresis model of the ECL file format.
This material law implements the hysteresis model of the ECL file format.
Definition: EclHysteresisTwoPhaseLaw.hpp:40
static void capillaryPressures(ContainerT &values, const Params ¶ms, const FluidState &fluidState)
Implements the multiplexer three phase capillary pressure law used by the ECLipse simulator...
Definition: EclMultiplexerMaterial.hpp:131
ParamsT Params
The type of the parameter objects for this law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:61
Implementation of a tabulated, piecewise linear capillary pressure law.
static void updateHysteresis(Params ¶ms, const FluidState &fluidState)
Update the hysteresis parameters after a time step.
Definition: EclMultiplexerMaterial.hpp:453
Scalar applySwatinit(unsigned elemIdx, Scalar pcow, Scalar Sw)
Modify the initial condition according to the SWATINIT keyword.
Definition: EclMaterialLawManager.hpp:156
Definition: Air_Mesitylene.hpp:33
const MaterialLawParams & connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const
Returns a material parameter object for a given element and saturation region.
Definition: EclMaterialLawManager.hpp:228
This file contains helper classes for the material laws.
Implementation of a tabulated, piecewise linear capillary pressure law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:52
This structure represents all values which can be possibly used as scaling points in the endpoint sca...
Definition: EclEpsScalingPoints.hpp:154
Implements a multiplexer class that provides all three phase capillary pressure laws used by the ECLi...
Definition: EclMultiplexerMaterial.hpp:57
This material law takes a material law defined for unscaled saturation and converts it to a material ...
Specifies the configuration used by the endpoint scaling code.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: SimpleModularFluidState.hpp:77
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
This material law takes a material law defined for unscaled saturation and converts it to a material ...
Definition: EclEpsTwoPhaseLaw.hpp:50
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.
Provides an simple way to create and manage the material law objects for a complete ECL deck...
Definition: EclMaterialLawManager.hpp:64
Specifies the configuration used by the ECL kr/pC hysteresis code.
A generic traits class for two-phase material laws.
Definition: MaterialTraits.hpp:59