28 #error "Eclipse input support in opm-common is required to use the ECL material manager!"
31 #ifndef OPM_ECL_MATERIAL_LAW_MANAGER_HPP
32 #define OPM_ECL_MATERIAL_LAW_MANAGER_HPP
46 #include <opm/common/OpmLog/OpmLog.hpp>
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>
67 template <
class TraitsT>
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 };
95 typedef typename GasOilEpsTwoPhaseLaw::Params GasOilEpsTwoPhaseParams;
96 typedef typename OilWaterEpsTwoPhaseLaw::Params OilWaterEpsTwoPhaseParams;
97 typedef typename GasWaterEpsTwoPhaseLaw::Params GasWaterEpsTwoPhaseParams;
103 typedef typename GasOilTwoPhaseLaw::Params GasOilTwoPhaseHystParams;
104 typedef typename OilWaterTwoPhaseLaw::Params OilWaterTwoPhaseHystParams;
105 typedef typename GasWaterTwoPhaseLaw::Params GasWaterTwoPhaseHystParams;
110 typedef typename MaterialLaw::Params MaterialLawParams;
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;
118 typedef std::vector<std::shared_ptr<EclEpsScalingPoints<Scalar> > > GasOilScalingPointsVector;
119 typedef std::vector<std::shared_ptr<EclEpsScalingPoints<Scalar> > > OilWaterScalingPointsVector;
120 typedef std::vector<std::shared_ptr<EclEpsScalingPoints<Scalar> > > GasWaterScalingPointsVector;
121 typedef std::vector<EclEpsScalingPointsInfo<Scalar>> OilWaterScalingInfoVector;
122 typedef std::vector<std::shared_ptr<GasOilTwoPhaseHystParams> > GasOilParamVector;
123 typedef std::vector<std::shared_ptr<OilWaterTwoPhaseHystParams> > OilWaterParamVector;
124 typedef std::vector<std::shared_ptr<GasWaterTwoPhaseHystParams> > GasWaterParamVector;
125 typedef std::vector<std::shared_ptr<MaterialLawParams> > MaterialLawParamsVector;
131 void initFromState(
const EclipseState& eclState)
134 const auto& runspec = eclState.runspec();
135 const size_t numSatRegions = runspec.tabdims().getNumSatTables();
137 const auto& ph = runspec.phases();
138 this->hasGas = ph.active(Phase::GAS);
139 this->hasOil = ph.active(Phase::OIL);
140 this->hasWater = ph.active(Phase::WATER);
142 readGlobalEpsOptions_(eclState);
143 readGlobalHysteresisOptions_(eclState);
144 readGlobalThreePhaseOptions_(runspec);
147 gasOilConfig = std::make_shared<EclEpsConfig>();
148 oilWaterConfig = std::make_shared<EclEpsConfig>();
149 gasWaterConfig = std::make_shared<EclEpsConfig>();
150 gasOilConfig->initFromState(eclState, EclGasOilSystem);
151 oilWaterConfig->initFromState(eclState, EclOilWaterSystem);
152 gasWaterConfig->initFromState(eclState, EclGasWaterSystem);
155 const auto& tables = eclState.getTableManager();
158 const auto& stone1exTables = tables.getStone1exTable();
160 if (! stone1exTables.empty()) {
162 stoneEtas.reserve(numSatRegions);
164 for (
const auto& table : stone1exTables) {
165 stoneEtas.push_back(table.eta);
170 this->unscaledEpsInfo_.resize(numSatRegions);
172 if (this->hasGas + this->hasOil + this->hasWater == 1) {
178 const auto tolcrit = runspec.saturationFunctionControls()
179 .minimumRelpermMobilityThreshold();
181 const auto rtep = satfunc::getRawTableEndpoints(tables, ph, tolcrit);
182 const auto rfunc = satfunc::getRawFunctionValues(tables, ph, rtep);
184 for (
unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) {
185 this->unscaledEpsInfo_[satRegionIdx]
186 .extractUnscaled(rtep, rfunc, satRegionIdx);
190 void initParamsForElements(
const EclipseState& eclState,
size_t numCompressedElems)
193 const size_t numSatRegions = eclState.runspec().tabdims().getNumSatTables();
196 gasOilUnscaledPointsVector_.resize(numSatRegions);
197 oilWaterUnscaledPointsVector_.resize(numSatRegions);
198 gasWaterUnscaledPointsVector_.resize(numSatRegions);
200 gasOilEffectiveParamVector_.resize(numSatRegions);
201 oilWaterEffectiveParamVector_.resize(numSatRegions);
202 gasWaterEffectiveParamVector_.resize(numSatRegions);
203 for (
unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) {
205 readGasOilUnscaledPoints_(gasOilUnscaledPointsVector_, gasOilConfig, eclState, satRegionIdx);
206 readOilWaterUnscaledPoints_(oilWaterUnscaledPointsVector_, oilWaterConfig, eclState, satRegionIdx);
207 readGasWaterUnscaledPoints_(gasWaterUnscaledPointsVector_, gasWaterConfig, eclState, satRegionIdx);
210 readGasOilEffectiveParameters_(gasOilEffectiveParamVector_, eclState, satRegionIdx);
211 readOilWaterEffectiveParameters_(oilWaterEffectiveParamVector_, eclState, satRegionIdx);
212 readGasWaterEffectiveParameters_(gasWaterEffectiveParamVector_, eclState, satRegionIdx);
217 satnumRegionArray_.resize(numCompressedElems);
218 if (eclState.fieldProps().has_int(
"SATNUM")) {
219 const auto& satnumRawData = eclState.fieldProps().get_int(
"SATNUM");
220 for (
unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
221 satnumRegionArray_[elemIdx] = satnumRawData[elemIdx] - 1;
225 std::fill(satnumRegionArray_.begin(), satnumRegionArray_.end(), 0);
229 imbnumRegionArray_ = satnumRegionArray_;
230 if (eclState.fieldProps().has_int(
"IMBNUM")) {
231 const auto& imbnumRawData = eclState.fieldProps().get_int(
"IMBNUM");
232 for (
unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
233 imbnumRegionArray_[elemIdx] = imbnumRawData[elemIdx] - 1;
237 assert(numCompressedElems == satnumRegionArray_.size());
238 assert(!enableHysteresis() || numCompressedElems == imbnumRegionArray_.size());
242 oilWaterScaledEpsInfoDrainage_.resize(numCompressedElems);
244 std::unique_ptr<EclEpsGridProperties> epsImbGridProperties;
246 if (enableHysteresis()) {
247 epsImbGridProperties = std::make_unique<EclEpsGridProperties>(eclState,
true);
251 materialLawParams_.resize(numCompressedElems);
253 for (
unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
254 unsigned satRegionIdx =
static_cast<unsigned>(satnumRegionArray_[elemIdx]);
255 auto gasOilParams = std::make_shared<GasOilTwoPhaseHystParams>();
256 auto oilWaterParams = std::make_shared<OilWaterTwoPhaseHystParams>();
257 auto gasWaterParams = std::make_shared<GasWaterTwoPhaseHystParams>();
258 gasOilParams->setConfig(hysteresisConfig_);
259 oilWaterParams->setConfig(hysteresisConfig_);
260 gasWaterParams->setConfig(hysteresisConfig_);
262 auto [gasOilScaledInfo, gasOilScaledPoint] =
263 readScaledPoints_(*gasOilConfig,
269 auto [owinfo, oilWaterScaledEpsPointDrainage] =
270 readScaledPoints_(*oilWaterConfig,
275 oilWaterScaledEpsInfoDrainage_[elemIdx] = owinfo;
277 auto [gasWaterScaledInfo, gasWaterScaledPoint] =
278 readScaledPoints_(*gasWaterConfig,
284 if (hasGas && hasOil) {
285 GasOilEpsTwoPhaseParams gasOilDrainParams;
286 gasOilDrainParams.setConfig(gasOilConfig);
287 gasOilDrainParams.setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
288 gasOilDrainParams.setScaledPoints(gasOilScaledPoint);
289 gasOilDrainParams.setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
290 gasOilDrainParams.finalize();
292 gasOilParams->setDrainageParams(gasOilDrainParams);
295 if (hasOil && hasWater) {
296 OilWaterEpsTwoPhaseParams oilWaterDrainParams;
297 oilWaterDrainParams.setConfig(oilWaterConfig);
298 oilWaterDrainParams.setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
299 oilWaterDrainParams.setScaledPoints(oilWaterScaledEpsPointDrainage);
300 oilWaterDrainParams.setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
301 oilWaterDrainParams.finalize();
303 oilWaterParams->setDrainageParams(oilWaterDrainParams);
306 if (hasGas && hasWater && !hasOil) {
307 GasWaterEpsTwoPhaseParams gasWaterDrainParams;
308 gasWaterDrainParams.setConfig(gasWaterConfig);
309 gasWaterDrainParams.setUnscaledPoints(gasWaterUnscaledPointsVector_[satRegionIdx]);
310 gasWaterDrainParams.setScaledPoints(gasWaterScaledPoint);
311 gasWaterDrainParams.setEffectiveLawParams(gasWaterEffectiveParamVector_[satRegionIdx]);
312 gasWaterDrainParams.finalize();
314 gasWaterParams->setDrainageParams(gasWaterDrainParams);
317 if (enableHysteresis()) {
318 auto [gasOilScaledImbInfo, gasOilScaledImbPoint] =
319 readScaledPoints_(*gasOilConfig,
321 *epsImbGridProperties,
325 auto [oilWaterScaledImbInfo, oilWaterScaledImbPoint] =
326 readScaledPoints_(*oilWaterConfig,
328 *epsImbGridProperties,
332 auto [gasWaterScaledImbInfo, gasWaterScaledImbPoint] =
333 readScaledPoints_(*gasWaterConfig,
335 *epsImbGridProperties,
339 unsigned imbRegionIdx = imbnumRegionArray_[elemIdx];
340 if (hasGas && hasOil) {
341 GasOilEpsTwoPhaseParams gasOilImbParamsHyst;
342 gasOilImbParamsHyst.setConfig(gasOilConfig);
343 gasOilImbParamsHyst.setUnscaledPoints(gasOilUnscaledPointsVector_[imbRegionIdx]);
344 gasOilImbParamsHyst.setScaledPoints(gasOilScaledImbPoint);
345 gasOilImbParamsHyst.setEffectiveLawParams(gasOilEffectiveParamVector_[imbRegionIdx]);
346 gasOilImbParamsHyst.finalize();
348 gasOilParams->setImbibitionParams(gasOilImbParamsHyst,
353 if (hasOil && hasWater) {
354 OilWaterEpsTwoPhaseParams oilWaterImbParamsHyst;
355 oilWaterImbParamsHyst.setConfig(oilWaterConfig);
356 oilWaterImbParamsHyst.setUnscaledPoints(oilWaterUnscaledPointsVector_[imbRegionIdx]);
357 oilWaterImbParamsHyst.setScaledPoints(oilWaterScaledImbPoint);
358 oilWaterImbParamsHyst.setEffectiveLawParams(oilWaterEffectiveParamVector_[imbRegionIdx]);
359 oilWaterImbParamsHyst.finalize();
361 oilWaterParams->setImbibitionParams(oilWaterImbParamsHyst,
366 if (hasGas && hasWater && !hasOil) {
367 GasWaterEpsTwoPhaseParams gasWaterImbParamsHyst;
368 gasWaterImbParamsHyst.setConfig(gasWaterConfig);
369 gasWaterImbParamsHyst.setUnscaledPoints(gasWaterUnscaledPointsVector_[imbRegionIdx]);
370 gasWaterImbParamsHyst.setScaledPoints(gasWaterScaledImbPoint);
371 gasWaterImbParamsHyst.setEffectiveLawParams(gasWaterEffectiveParamVector_[imbRegionIdx]);
372 gasWaterImbParamsHyst.finalize();
374 gasWaterParams->setImbibitionParams(gasWaterImbParamsHyst,
375 gasWaterScaledImbInfo,
380 if (hasGas && hasOil)
381 gasOilParams->finalize();
383 if (hasOil && hasWater)
384 oilWaterParams->finalize();
386 if (hasGas && hasWater && !hasOil)
387 gasWaterParams->finalize();
389 initThreePhaseParams_(eclState,
390 materialLawParams_[elemIdx],
392 oilWaterScaledEpsInfoDrainage_[elemIdx],
397 materialLawParams_[elemIdx].finalize();
414 auto& elemScaledEpsInfo = oilWaterScaledEpsInfoDrainage_[elemIdx];
419 Sw = elemScaledEpsInfo.Swu;
422 if (Sw <= elemScaledEpsInfo.Swl)
423 Sw = elemScaledEpsInfo.Swl;
439 fs.setSaturation(waterPhaseIdx, Sw);
440 fs.setSaturation(gasPhaseIdx, 0);
441 fs.setSaturation(oilPhaseIdx, 0);
442 Scalar pc[numPhases] = { 0 };
445 Scalar pcowAtSw = pc[oilPhaseIdx] - pc[waterPhaseIdx];
446 const Scalar pcowAtSwThreshold = 1.0;
448 if (std::abs(pcowAtSw) > pcowAtSwThreshold) {
449 elemScaledEpsInfo.maxPcow *= pcow/pcowAtSw;
450 auto& elemEclEpsScalingPoints = oilWaterScaledEpsPointsDrainage(elemIdx);
451 elemEclEpsScalingPoints.init(elemScaledEpsInfo, *oilWaterEclEpsConfig_, EclOilWaterSystem);
458 bool enableEndPointScaling()
const
459 {
return enableEndPointScaling_; }
461 bool enableHysteresis()
const
462 {
return hysteresisConfig_->enableHysteresis(); }
464 MaterialLawParams& materialLawParams(
unsigned elemIdx)
466 assert(elemIdx < materialLawParams_.size());
467 return materialLawParams_[elemIdx];
470 const MaterialLawParams& materialLawParams(
unsigned elemIdx)
const
472 assert(elemIdx < materialLawParams_.size());
473 return materialLawParams_[elemIdx];
486 MaterialLawParams& mlp =
const_cast<MaterialLawParams&
>(materialLawParams_[elemIdx]);
489 if (enableHysteresis())
490 OpmLog::warning(
"Warning: Using non-default satnum regions for connection is not tested in combination with hysteresis");
496 switch (mlp.approach()) {
497 case EclMultiplexerApproach::EclStone1Approach: {
498 auto& realParams = mlp.template getRealParams<EclMultiplexerApproach::EclStone1Approach>();
500 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
501 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
502 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
503 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
513 case EclMultiplexerApproach::EclStone2Approach: {
514 auto& realParams = mlp.template getRealParams<EclMultiplexerApproach::EclStone2Approach>();
515 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
516 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
517 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
518 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
528 case EclMultiplexerApproach::EclDefaultApproach: {
529 auto& realParams = mlp.template getRealParams<EclMultiplexerApproach::EclDefaultApproach>();
530 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
531 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
532 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
533 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
543 case EclMultiplexerApproach::EclTwoPhaseApproach: {
544 auto& realParams = mlp.template getRealParams<EclMultiplexerApproach::EclTwoPhaseApproach>();
545 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
546 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
547 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
548 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
559 throw std::logic_error(
"Enum value for material approach unknown!");
565 int satnumRegionIdx(
unsigned elemIdx)
const
566 {
return satnumRegionArray_[elemIdx]; }
568 int imbnumRegionIdx(
unsigned elemIdx)
const
569 {
return imbnumRegionArray_[elemIdx]; }
571 std::shared_ptr<MaterialLawParams>& materialLawParamsPointerReferenceHack(
unsigned elemIdx)
573 assert(0 <= elemIdx && elemIdx < materialLawParams_.size());
574 return materialLawParams_[elemIdx];
577 template <
class Flu
idState>
578 void updateHysteresis(
const FluidState& fluidState,
unsigned elemIdx)
580 if (!enableHysteresis())
586 void oilWaterHysteresisParams(Scalar& pcSwMdc,
588 unsigned elemIdx)
const
590 if (!enableHysteresis())
591 throw std::runtime_error(
"Cannot get hysteresis parameters if hysteresis not enabled.");
593 const auto& params = materialLawParams(elemIdx);
594 MaterialLaw::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, params);
597 void setOilWaterHysteresisParams(
const Scalar& pcSwMdc,
598 const Scalar& krnSwMdc,
601 if (!enableHysteresis())
602 throw std::runtime_error(
"Cannot set hysteresis parameters if hysteresis not enabled.");
604 auto& params = materialLawParams(elemIdx);
605 MaterialLaw::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, params);
608 void gasOilHysteresisParams(Scalar& pcSwMdc,
610 unsigned elemIdx)
const
612 if (!enableHysteresis())
613 throw std::runtime_error(
"Cannot get hysteresis parameters if hysteresis not enabled.");
615 const auto& params = materialLawParams(elemIdx);
616 MaterialLaw::gasOilHysteresisParams(pcSwMdc, krnSwMdc, params);
619 void setGasOilHysteresisParams(
const Scalar& pcSwMdc,
620 const Scalar& krnSwMdc,
623 if (!enableHysteresis())
624 throw std::runtime_error(
"Cannot set hysteresis parameters if hysteresis not enabled.");
626 auto& params = materialLawParams(elemIdx);
627 MaterialLaw::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, params);
630 EclEpsScalingPoints<Scalar>& oilWaterScaledEpsPointsDrainage(
unsigned elemIdx)
632 auto& materialParams = materialLawParams_[elemIdx];
633 switch (materialParams.approach()) {
634 case EclMultiplexerApproach::EclStone1Approach: {
635 auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclStone1Approach>();
636 return realParams.oilWaterParams().drainageParams().scaledPoints();
639 case EclMultiplexerApproach::EclStone2Approach: {
640 auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclStone2Approach>();
641 return realParams.oilWaterParams().drainageParams().scaledPoints();
644 case EclMultiplexerApproach::EclDefaultApproach: {
645 auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclDefaultApproach>();
646 return realParams.oilWaterParams().drainageParams().scaledPoints();
649 case EclMultiplexerApproach::EclTwoPhaseApproach: {
650 auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclTwoPhaseApproach>();
651 return realParams.oilWaterParams().drainageParams().scaledPoints();
654 throw std::logic_error(
"Enum value for material approach unknown!");
658 const EclEpsScalingPointsInfo<Scalar>& oilWaterScaledEpsInfoDrainage(
size_t elemIdx)
const
659 {
return oilWaterScaledEpsInfoDrainage_[elemIdx]; }
662 void readGlobalEpsOptions_(
const EclipseState& eclState)
664 oilWaterEclEpsConfig_ = std::make_shared<EclEpsConfig>();
665 oilWaterEclEpsConfig_->initFromState(eclState, EclOilWaterSystem);
667 enableEndPointScaling_ = eclState.getTableManager().hasTables(
"ENKRVD");
670 void readGlobalHysteresisOptions_(
const EclipseState& state)
672 hysteresisConfig_ = std::make_shared<EclHysteresisConfig>();
673 hysteresisConfig_->initFromState(state.runspec());
676 void readGlobalThreePhaseOptions_(
const Runspec& runspec)
678 bool gasEnabled = runspec.phases().active(Phase::GAS);
679 bool oilEnabled = runspec.phases().active(Phase::OIL);
680 bool waterEnabled = runspec.phases().active(Phase::WATER);
685 + (waterEnabled?1:0);
687 if (numEnabled == 0) {
688 throw std::runtime_error(
"At least one fluid phase must be enabled. (Is: "+std::to_string(numEnabled)+
")");
689 }
else if (numEnabled == 1) {
690 threePhaseApproach_ = EclMultiplexerApproach::EclOnePhaseApproach;
691 }
else if ( numEnabled == 2) {
692 threePhaseApproach_ = EclMultiplexerApproach::EclTwoPhaseApproach;
694 twoPhaseApproach_ = EclTwoPhaseApproach::EclTwoPhaseOilWater;
695 else if (!oilEnabled)
696 twoPhaseApproach_ = EclTwoPhaseApproach::EclTwoPhaseGasWater;
697 else if (!waterEnabled)
698 twoPhaseApproach_ = EclTwoPhaseApproach::EclTwoPhaseGasOil;
701 assert(numEnabled == 3);
703 threePhaseApproach_ = EclMultiplexerApproach::EclDefaultApproach;
704 const auto& satctrls = runspec.saturationFunctionControls();
705 if (satctrls.krModel() == SatFuncControls::ThreePhaseOilKrModel::Stone2)
706 threePhaseApproach_ = EclMultiplexerApproach::EclStone2Approach;
707 else if (satctrls.krModel() == SatFuncControls::ThreePhaseOilKrModel::Stone1)
708 threePhaseApproach_ = EclMultiplexerApproach::EclStone1Approach;
712 template <
class Container>
713 void readGasOilEffectiveParameters_(Container& dest,
714 const EclipseState& eclState,
715 unsigned satRegionIdx)
717 if (!hasGas || !hasOil)
721 dest[satRegionIdx] = std::make_shared<GasOilEffectiveTwoPhaseParams>();
723 auto& effParams = *dest[satRegionIdx];
727 const Scalar Swco = unscaledEpsInfo_[satRegionIdx].Swl;
728 const auto tolcrit = eclState.runspec().saturationFunctionControls()
729 .minimumRelpermMobilityThreshold();
731 const auto& tableManager = eclState.getTableManager();
733 switch (eclState.runspec().saturationFunctionControls().family()) {
734 case SatFuncControls::KeywordFamily::Family_I:
736 const TableContainer& sgofTables = tableManager.getSgofTables();
737 const TableContainer& slgofTables = tableManager.getSlgofTables();
738 if (!sgofTables.empty())
739 readGasOilEffectiveParametersSgof_(effParams, Swco, tolcrit,
740 sgofTables.getTable<SgofTable>(satRegionIdx));
741 else if (!slgofTables.empty())
742 readGasOilEffectiveParametersSlgof_(effParams, Swco, tolcrit,
743 slgofTables.getTable<SlgofTable>(satRegionIdx));
747 case SatFuncControls::KeywordFamily::Family_II:
749 const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable<SgfnTable>( satRegionIdx );
752 const Sof2Table& sof2Table = tableManager.getSof2Tables().getTable<Sof2Table>( satRegionIdx );
753 readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof2Table, sgfnTable);
756 const Sof3Table& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>( satRegionIdx );
757 readGasOilEffectiveParametersFamily2_(effParams, Swco, tolcrit, sof3Table, sgfnTable);
762 case SatFuncControls::KeywordFamily::Undefined:
763 throw std::domain_error(
"No valid saturation keyword family specified");
767 void readGasOilEffectiveParametersSgof_(GasOilEffectiveTwoPhaseParams& effParams,
769 const double tolcrit,
770 const SgofTable& sgofTable)
773 std::vector<double> SoSamples(sgofTable.numRows());
774 for (
size_t sampleIdx = 0; sampleIdx < sgofTable.numRows(); ++ sampleIdx) {
775 SoSamples[sampleIdx] = (1.0 - Swco) - sgofTable.get(
"SG", sampleIdx);
778 effParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn(
"KROG")));
779 effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgofTable.getColumn(
"KRG")));
780 effParams.setPcnwSamples(SoSamples, sgofTable.getColumn(
"PCOG").vectorCopy());
781 effParams.finalize();
784 void readGasOilEffectiveParametersSlgof_(GasOilEffectiveTwoPhaseParams& effParams,
786 const double tolcrit,
787 const SlgofTable& slgofTable)
790 std::vector<double> SoSamples(slgofTable.numRows());
791 for (
size_t sampleIdx = 0; sampleIdx < slgofTable.numRows(); ++ sampleIdx) {
792 SoSamples[sampleIdx] = slgofTable.get(
"SL", sampleIdx) - Swco;
795 effParams.setKrwSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn(
"KROG")));
796 effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, slgofTable.getColumn(
"KRG")));
797 effParams.setPcnwSamples(SoSamples, slgofTable.getColumn(
"PCOG").vectorCopy());
798 effParams.finalize();
801 void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
803 const double tolcrit,
804 const Sof3Table& sof3Table,
805 const SgfnTable& sgfnTable)
808 std::vector<double> SoSamples(sgfnTable.numRows());
809 std::vector<double> SoColumn = sof3Table.getColumn(
"SO").vectorCopy();
810 for (
size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) {
811 SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get(
"SG", sampleIdx);
814 effParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof3Table.getColumn(
"KROG")));
815 effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn(
"KRG")));
816 effParams.setPcnwSamples(SoSamples, sgfnTable.getColumn(
"PCOG").vectorCopy());
817 effParams.finalize();
820 void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
822 const double tolcrit,
823 const Sof2Table& sof2Table,
824 const SgfnTable& sgfnTable)
827 std::vector<double> SoSamples(sgfnTable.numRows());
828 std::vector<double> SoColumn = sof2Table.getColumn(
"SO").vectorCopy();
829 for (
size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) {
830 SoSamples[sampleIdx] = (1.0 - Swco) - sgfnTable.get(
"SG", sampleIdx);
833 effParams.setKrwSamples(SoColumn, normalizeKrValues_(tolcrit, sof2Table.getColumn(
"KRO")));
834 effParams.setKrnSamples(SoSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn(
"KRG")));
835 effParams.setPcnwSamples(SoSamples, sgfnTable.getColumn(
"PCOG").vectorCopy());
836 effParams.finalize();
839 template <
class Container>
840 void readOilWaterEffectiveParameters_(Container& dest,
841 const EclipseState& eclState,
842 unsigned satRegionIdx)
844 if (!hasOil || !hasWater)
848 dest[satRegionIdx] = std::make_shared<OilWaterEffectiveTwoPhaseParams>();
850 const auto tolcrit = eclState.runspec().saturationFunctionControls()
851 .minimumRelpermMobilityThreshold();
853 const auto& tableManager = eclState.getTableManager();
854 auto& effParams = *dest[satRegionIdx];
856 switch (eclState.runspec().saturationFunctionControls().family()) {
857 case SatFuncControls::KeywordFamily::Family_I:
859 const auto& swofTable = tableManager.getSwofTables().getTable<SwofTable>(satRegionIdx);
860 const std::vector<double> SwColumn = swofTable.getColumn(
"SW").vectorCopy();
862 effParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn(
"KRW")));
863 effParams.setKrnSamples(SwColumn, normalizeKrValues_(tolcrit, swofTable.getColumn(
"KROW")));
864 effParams.setPcnwSamples(SwColumn, swofTable.getColumn(
"PCOW").vectorCopy());
865 effParams.finalize();
869 case SatFuncControls::KeywordFamily::Family_II:
871 const auto& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>(satRegionIdx);
872 const std::vector<double> SwColumn = swfnTable.getColumn(
"SW").vectorCopy();
873 effParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn(
"KRW")));
874 effParams.setPcnwSamples(SwColumn, swfnTable.getColumn(
"PCOW").vectorCopy());
877 const auto& sof2Table = tableManager.getSof2Tables().getTable<Sof2Table>(satRegionIdx);
879 std::vector<double> SwSamples(sof2Table.numRows());
880 for (
size_t sampleIdx = 0; sampleIdx < sof2Table.numRows(); ++ sampleIdx)
881 SwSamples[sampleIdx] = 1 - sof2Table.get(
"SO", sampleIdx);
883 effParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof2Table.getColumn(
"KRO")));
885 const auto& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>(satRegionIdx);
887 std::vector<double> SwSamples(sof3Table.numRows());
888 for (
size_t sampleIdx = 0; sampleIdx < sof3Table.numRows(); ++ sampleIdx)
889 SwSamples[sampleIdx] = 1 - sof3Table.get(
"SO", sampleIdx);
891 effParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sof3Table.getColumn(
"KROW")));
893 effParams.finalize();
897 case SatFuncControls::KeywordFamily::Undefined:
898 throw std::domain_error(
"No valid saturation keyword family specified");
902 template <
class Container>
903 void readGasWaterEffectiveParameters_(Container& dest,
904 const EclipseState& eclState,
905 unsigned satRegionIdx)
907 if (!hasGas || !hasWater || hasOil)
911 dest[satRegionIdx] = std::make_shared<GasWaterEffectiveTwoPhaseParams>();
913 auto& effParams = *dest[satRegionIdx];
915 const auto tolcrit = eclState.runspec().saturationFunctionControls()
916 .minimumRelpermMobilityThreshold();
918 const auto& tableManager = eclState.getTableManager();
920 switch (eclState.runspec().saturationFunctionControls().family()) {
921 case SatFuncControls::KeywordFamily::Family_I:
923 throw std::domain_error(
"Saturation keyword family I is not applicable for a gas-water system");
926 case SatFuncControls::KeywordFamily::Family_II:
929 const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable<SgfnTable>( satRegionIdx );
930 const SwfnTable& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>( satRegionIdx );
932 std::vector<double> SwColumn = swfnTable.getColumn(
"SW").vectorCopy();
934 effParams.setKrwSamples(SwColumn, normalizeKrValues_(tolcrit, swfnTable.getColumn(
"KRW")));
935 std::vector<double> SwSamples(sgfnTable.numRows());
936 for (
size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx)
937 SwSamples[sampleIdx] = 1 - sgfnTable.get(
"SG", sampleIdx);
938 effParams.setKrnSamples(SwSamples, normalizeKrValues_(tolcrit, sgfnTable.getColumn(
"KRG")));
941 effParams.setPcnwSamples(SwColumn, swfnTable.getColumn(
"PCOW").vectorCopy());
942 effParams.finalize();
947 case SatFuncControls::KeywordFamily::Undefined:
948 throw std::domain_error(
"No valid saturation keyword family specified");
952 template <
class Container>
953 void readGasOilUnscaledPoints_(Container& dest,
954 std::shared_ptr<EclEpsConfig> config,
955 const EclipseState& ,
956 unsigned satRegionIdx)
958 if (!hasGas || !hasOil)
962 dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
963 dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], *config, EclGasOilSystem);
966 template <
class Container>
967 void readOilWaterUnscaledPoints_(Container& dest,
968 std::shared_ptr<EclEpsConfig> config,
969 const EclipseState& ,
970 unsigned satRegionIdx)
972 if (!hasOil || !hasWater)
976 dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
977 dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], *config, EclOilWaterSystem);
980 template <
class Container>
981 void readGasWaterUnscaledPoints_(Container& dest,
982 std::shared_ptr<EclEpsConfig> config,
983 const EclipseState& ,
984 unsigned satRegionIdx)
990 dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
991 dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], *config, EclGasWaterSystem);
994 std::tuple<EclEpsScalingPointsInfo<Scalar>,
995 EclEpsScalingPoints<Scalar>>
996 readScaledPoints_(
const EclEpsConfig& config,
997 const EclipseState& eclState,
998 const EclEpsGridProperties& epsGridProperties,
1000 EclTwoPhaseSystemType type)
1002 unsigned satRegionIdx = epsGridProperties.satRegion( elemIdx );
1004 EclEpsScalingPointsInfo<Scalar> destInfo(unscaledEpsInfo_[satRegionIdx]);
1005 destInfo.extractScaled(eclState, epsGridProperties, elemIdx);
1007 EclEpsScalingPoints<Scalar> destPoint;
1008 destPoint.init(destInfo, config, type);
1010 return {destInfo, destPoint};
1013 void initThreePhaseParams_(
const EclipseState& ,
1014 MaterialLawParams& materialParams,
1015 unsigned satRegionIdx,
1016 const EclEpsScalingPointsInfo<Scalar>& epsInfo,
1017 std::shared_ptr<OilWaterTwoPhaseHystParams> oilWaterParams,
1018 std::shared_ptr<GasOilTwoPhaseHystParams> gasOilParams,
1019 std::shared_ptr<GasWaterTwoPhaseHystParams> gasWaterParams)
1021 materialParams.setApproach(threePhaseApproach_);
1023 switch (materialParams.approach()) {
1024 case EclMultiplexerApproach::EclStone1Approach: {
1025 auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclStone1Approach>();
1026 realParams.setGasOilParams(gasOilParams);
1027 realParams.setOilWaterParams(oilWaterParams);
1028 realParams.setSwl(epsInfo.Swl);
1030 if (!stoneEtas.empty()) {
1031 realParams.setEta(stoneEtas[satRegionIdx]);
1034 realParams.setEta(1.0);
1035 realParams.finalize();
1039 case EclMultiplexerApproach::EclStone2Approach: {
1040 auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclStone2Approach>();
1041 realParams.setGasOilParams(gasOilParams);
1042 realParams.setOilWaterParams(oilWaterParams);
1043 realParams.setSwl(epsInfo.Swl);
1044 realParams.finalize();
1048 case EclMultiplexerApproach::EclDefaultApproach: {
1049 auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclDefaultApproach>();
1050 realParams.setGasOilParams(gasOilParams);
1051 realParams.setOilWaterParams(oilWaterParams);
1052 realParams.setSwl(epsInfo.Swl);
1053 realParams.finalize();
1057 case EclMultiplexerApproach::EclTwoPhaseApproach: {
1058 auto& realParams = materialParams.template getRealParams<EclMultiplexerApproach::EclTwoPhaseApproach>();
1059 realParams.setGasOilParams(gasOilParams);
1060 realParams.setOilWaterParams(oilWaterParams);
1061 realParams.setGasWaterParams(gasWaterParams);
1062 realParams.setApproach(twoPhaseApproach_);
1063 realParams.finalize();
1067 case EclMultiplexerApproach::EclOnePhaseApproach: {
1075 std::vector<double> normalizeKrValues_(
const double tolcrit,
1076 const TableColumn& krValues)
const
1078 auto kr = krValues.vectorCopy();
1079 std::transform(kr.begin(), kr.end(), kr.begin(),
1080 [tolcrit](
const double kri)
1082 return (kri > tolcrit) ? kri : 0.0;
1088 bool enableEndPointScaling_;
1089 std::shared_ptr<EclHysteresisConfig> hysteresisConfig_;
1091 std::shared_ptr<EclEpsConfig> oilWaterEclEpsConfig_;
1092 std::vector<EclEpsScalingPointsInfo<Scalar>> unscaledEpsInfo_;
1093 OilWaterScalingInfoVector oilWaterScaledEpsInfoDrainage_;
1095 std::shared_ptr<EclEpsConfig> gasWaterEclEpsConfig_;
1097 GasOilScalingPointsVector gasOilUnscaledPointsVector_;
1098 OilWaterScalingPointsVector oilWaterUnscaledPointsVector_;
1099 GasWaterScalingPointsVector gasWaterUnscaledPointsVector_;
1101 GasOilEffectiveParamVector gasOilEffectiveParamVector_;
1102 OilWaterEffectiveParamVector oilWaterEffectiveParamVector_;
1103 GasWaterEffectiveParamVector gasWaterEffectiveParamVector_;
1105 EclMultiplexerApproach threePhaseApproach_ = EclMultiplexerApproach::EclDefaultApproach;
1107 enum EclTwoPhaseApproach twoPhaseApproach_ = EclTwoPhaseApproach::EclTwoPhaseGasOil;
1109 std::vector<MaterialLawParams> materialLawParams_;
1111 std::vector<int> satnumRegionArray_;
1112 std::vector<int> imbnumRegionArray_;
1113 std::vector<Scalar> stoneEtas;
1119 std::shared_ptr<EclEpsConfig> gasOilConfig;
1120 std::shared_ptr<EclEpsConfig> oilWaterConfig;
1121 std::shared_ptr<EclEpsConfig> gasWaterConfig;
Specifies the configuration used by the endpoint scaling code.
This material law takes a material law defined for unscaled saturation and converts it to a material ...
Specifies the configuration used by the ECL kr/pC hysteresis code.
This material law implements the hysteresis model of the ECL file format.
Implements a multiplexer class that provides all three phase capillary pressure laws used by the ECLi...
Implementation for the parameters required by the material law for two-phase simulations.
This file contains helper classes for the material laws.
Implementation of a tabulated, piecewise linear capillary pressure law.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: EclEpsGridProperties.hpp:76
This material law takes a material law defined for unscaled saturation and converts it to a material ...
Definition: EclEpsTwoPhaseLaw.hpp:53
This material law implements the hysteresis model of the ECL file format.
Definition: EclHysteresisTwoPhaseLaw.hpp:41
Provides an simple way to create and manage the material law objects for a complete ECL deck.
Definition: EclMaterialLawManager.hpp:69
const MaterialLawParams & connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const
Returns a material parameter object for a given element and saturation region.
Definition: EclMaterialLawManager.hpp:484
Scalar applySwatinit(unsigned elemIdx, Scalar pcow, Scalar Sw)
Modify the initial condition according to the SWATINIT keyword.
Definition: EclMaterialLawManager.hpp:410
Implements a multiplexer class that provides all three phase capillary pressure laws used by the ECLi...
Definition: EclMultiplexerMaterial.hpp:58
static void capillaryPressures(ContainerT &values, const Params ¶ms, const FluidState &fluidState)
Implements the multiplexer three phase capillary pressure law used by the ECLipse simulator.
Definition: EclMultiplexerMaterial.hpp:135
static void updateHysteresis(Params ¶ms, 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