00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00027 #if ! HAVE_OPM_PARSER
00028 #error "The opm-parser module is required to use the ECL material manager!"
00029 #endif
00030
00031 #ifndef OPM_ECL_MATERIAL_LAW_MANAGER_HPP
00032 #define OPM_ECL_MATERIAL_LAW_MANAGER_HPP
00033
00034 #include <opm/material/fluidmatrixinteractions/EclTwoPhaseMaterialParams.hpp>
00035 #include <opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp>
00036 #include <opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp>
00037 #include <opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLaw.hpp>
00038 #include <opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp>
00039 #include <opm/material/fluidmatrixinteractions/EclEpsConfig.hpp>
00040 #include <opm/material/fluidmatrixinteractions/EclHysteresisConfig.hpp>
00041 #include <opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp>
00042 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
00043 #include <opm/material/fluidstates/SimpleModularFluidState.hpp>
00044
00045 #include <opm/common/Exceptions.hpp>
00046 #include <opm/common/ErrorMacros.hpp>
00047
00048 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
00049 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
00050 #include <opm/parser/eclipse/Deck/Deck.hpp>
00051
00052 #include <algorithm>
00053
00054
00055 namespace Opm {
00056
00063 template <class TraitsT>
00064 class EclMaterialLawManager
00065 {
00066 private:
00067 typedef TraitsT Traits;
00068 typedef typename Traits::Scalar Scalar;
00069 enum { waterPhaseIdx = Traits::wettingPhaseIdx };
00070 enum { oilPhaseIdx = Traits::nonWettingPhaseIdx };
00071 enum { gasPhaseIdx = Traits::gasPhaseIdx };
00072 enum { numPhases = Traits::numPhases };
00073
00074 typedef TwoPhaseMaterialTraits<Scalar, oilPhaseIdx, gasPhaseIdx> GasOilTraits;
00075 typedef TwoPhaseMaterialTraits<Scalar, waterPhaseIdx, oilPhaseIdx> OilWaterTraits;
00076
00077
00078 typedef PiecewiseLinearTwoPhaseMaterial<GasOilTraits> GasOilEffectiveTwoPhaseLaw;
00079 typedef PiecewiseLinearTwoPhaseMaterial<OilWaterTraits> OilWaterEffectiveTwoPhaseLaw;
00080 typedef typename GasOilEffectiveTwoPhaseLaw::Params GasOilEffectiveTwoPhaseParams;
00081 typedef typename OilWaterEffectiveTwoPhaseLaw::Params OilWaterEffectiveTwoPhaseParams;
00082
00083
00084 typedef EclEpsTwoPhaseLaw<GasOilEffectiveTwoPhaseLaw> GasOilEpsTwoPhaseLaw;
00085 typedef EclEpsTwoPhaseLaw<OilWaterEffectiveTwoPhaseLaw> OilWaterEpsTwoPhaseLaw;
00086 typedef typename GasOilEpsTwoPhaseLaw::Params GasOilEpsTwoPhaseParams;
00087 typedef typename OilWaterEpsTwoPhaseLaw::Params OilWaterEpsTwoPhaseParams;
00088
00089
00090 typedef EclHysteresisTwoPhaseLaw<GasOilEpsTwoPhaseLaw> GasOilTwoPhaseLaw;
00091 typedef EclHysteresisTwoPhaseLaw<OilWaterEpsTwoPhaseLaw> OilWaterTwoPhaseLaw;
00092 typedef typename GasOilTwoPhaseLaw::Params GasOilTwoPhaseHystParams;
00093 typedef typename OilWaterTwoPhaseLaw::Params OilWaterTwoPhaseHystParams;
00094
00095 public:
00096
00097 typedef EclMultiplexerMaterial<Traits, GasOilTwoPhaseLaw, OilWaterTwoPhaseLaw> MaterialLaw;
00098 typedef typename MaterialLaw::Params MaterialLawParams;
00099
00100 private:
00101
00102 typedef std::vector<std::shared_ptr<GasOilEffectiveTwoPhaseParams> > GasOilEffectiveParamVector;
00103 typedef std::vector<std::shared_ptr<OilWaterEffectiveTwoPhaseParams> > OilWaterEffectiveParamVector;
00104 typedef std::vector<std::shared_ptr<EclEpsScalingPoints<Scalar> > > GasOilScalingPointsVector;
00105 typedef std::vector<std::shared_ptr<EclEpsScalingPoints<Scalar> > > OilWaterScalingPointsVector;
00106 typedef std::vector<std::shared_ptr<EclEpsScalingPointsInfo<Scalar> > > GasOilScalingInfoVector;
00107 typedef std::vector<std::shared_ptr<EclEpsScalingPointsInfo<Scalar> > > OilWaterScalingInfoVector;
00108 typedef std::vector<std::shared_ptr<GasOilTwoPhaseHystParams> > GasOilParamVector;
00109 typedef std::vector<std::shared_ptr<OilWaterTwoPhaseHystParams> > OilWaterParamVector;
00110 typedef std::vector<std::shared_ptr<MaterialLawParams> > MaterialLawParamsVector;
00111
00112 public:
00113 EclMaterialLawManager()
00114 {}
00115
00116 void initFromDeck(const Opm::Deck& deck,
00117 const Opm::EclipseState& eclState,
00118 const std::vector<int>& compressedToCartesianElemIdx)
00119 {
00120
00121 const size_t numSatRegions = eclState.runspec().tabdims().getNumSatTables();
00122 size_t numCompressedElems = compressedToCartesianElemIdx.size();
00123
00124
00125
00126 satnumRegionArray_.resize(numCompressedElems);
00127 if (eclState.get3DProperties().hasDeckIntGridProperty("SATNUM")) {
00128 const auto& satnumRawData = eclState.get3DProperties().getIntGridProperty("SATNUM").getData();
00129 for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
00130 unsigned cartesianElemIdx = static_cast<unsigned>(compressedToCartesianElemIdx[elemIdx]);
00131 satnumRegionArray_[elemIdx] = satnumRawData[cartesianElemIdx] - 1;
00132 }
00133 }
00134 else
00135 std::fill(satnumRegionArray_.begin(), satnumRegionArray_.end(), 0);
00136
00137 readGlobalEpsOptions_(deck, eclState);
00138 readGlobalHysteresisOptions_(deck);
00139 readGlobalThreePhaseOptions_(deck);
00140
00141 unscaledEpsInfo_.resize(numSatRegions);
00142 for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx)
00143 unscaledEpsInfo_[satRegionIdx].extractUnscaled(deck, eclState, satRegionIdx);
00144
00145 initParamsForElements_(deck, eclState, compressedToCartesianElemIdx, satnumRegionArray_);
00146 }
00147
00156 Scalar applySwatinit(unsigned elemIdx,
00157 Scalar pcow,
00158 Scalar Sw)
00159 {
00160 auto& elemScaledEpsInfo = *oilWaterScaledEpsInfoDrainage_[elemIdx];
00161
00162
00163
00164 if (pcow < 0.0)
00165 Sw = elemScaledEpsInfo.Swu;
00166 else {
00167
00168 if (Sw <= elemScaledEpsInfo.Swl)
00169 Sw = elemScaledEpsInfo.Swl;
00170
00171
00172 typedef Opm::SimpleModularFluidState<Scalar,
00173 numPhases,
00174 0,
00175 void,
00176 false,
00177 false,
00178 false,
00179 false,
00180 true,
00181 false,
00182 false,
00183 false> FluidState;
00184 FluidState fs;
00185 fs.setSaturation(waterPhaseIdx, Sw);
00186 fs.setSaturation(gasPhaseIdx, 0);
00187 fs.setSaturation(oilPhaseIdx, 0);
00188 Scalar pc[numPhases] = { 0 };
00189 MaterialLaw::capillaryPressures(pc, materialLawParams(elemIdx), fs);
00190
00191 Scalar pcowAtSw = pc[oilPhaseIdx] - pc[waterPhaseIdx];
00192 if (pcowAtSw > 0.0) {
00193 elemScaledEpsInfo.maxPcow *= pcow/pcowAtSw;
00194 auto& elemEclEpsScalingPoints = oilWaterScaledEpsPointsDrainage(elemIdx);
00195 elemEclEpsScalingPoints.init(elemScaledEpsInfo, *oilWaterEclEpsConfig_, Opm::EclOilWaterSystem);
00196 }
00197 }
00198
00199 return Sw;
00200 }
00201
00202 bool enableEndPointScaling() const
00203 { return enableEndPointScaling_; }
00204
00205 bool enableHysteresis() const
00206 { return hysteresisConfig_->enableHysteresis(); }
00207
00208 MaterialLawParams& materialLawParams(unsigned elemIdx)
00209 {
00210 assert(0 <= elemIdx && elemIdx < materialLawParams_.size());
00211 return *materialLawParams_[elemIdx];
00212 }
00213
00214 const MaterialLawParams& materialLawParams(unsigned elemIdx) const
00215 {
00216 assert(0 <= elemIdx && elemIdx < materialLawParams_.size());
00217 return *materialLawParams_[elemIdx];
00218 }
00219
00228 const MaterialLawParams& connectionMaterialLawParams(unsigned satRegionIdx, unsigned elemIdx) const
00229 {
00230 MaterialLawParams& mlp = *materialLawParams_[elemIdx];
00231
00232 if (enableHysteresis()) {
00233 OPM_MESSAGE("Warning: Using non-defautl satnum regions for conenction is not tested in combination with hysteresis");
00234 }
00235
00236
00237
00238
00239 switch (mlp.approach()) {
00240 case EclStone1Approach: {
00241 auto& realParams = mlp.template getRealParams<Opm::EclStone1Approach>();
00242
00243 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
00244 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
00245 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
00246 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
00247
00248
00249
00250
00251
00252
00253 }
00254 break;
00255
00256 case EclStone2Approach: {
00257 auto& realParams = mlp.template getRealParams<Opm::EclStone2Approach>();
00258 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
00259 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
00260 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
00261 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
00262
00263
00264
00265
00266
00267
00268 }
00269 break;
00270
00271 case EclDefaultApproach: {
00272 auto& realParams = mlp.template getRealParams<Opm::EclDefaultApproach>();
00273 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
00274 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
00275 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
00276 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
00277
00278
00279
00280
00281
00282
00283 }
00284 break;
00285
00286 case EclTwoPhaseApproach: {
00287 auto& realParams = mlp.template getRealParams<Opm::EclTwoPhaseApproach>();
00288 realParams.oilWaterParams().drainageParams().setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
00289 realParams.oilWaterParams().drainageParams().setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
00290 realParams.gasOilParams().drainageParams().setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
00291 realParams.gasOilParams().drainageParams().setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
00292
00293
00294
00295
00296
00297
00298 }
00299 break;
00300
00301 default:
00302 OPM_THROW(std::logic_error, "Enum value for material approach unknown!");
00303 }
00304
00305 return mlp;
00306 }
00307
00308 int satnumRegionIdx(unsigned elemIdx) const {
00309 return satnumRegionArray_[elemIdx];
00310 }
00311
00312 std::shared_ptr<MaterialLawParams>& materialLawParamsPointerReferenceHack(unsigned elemIdx)
00313 {
00314 assert(0 <= elemIdx && elemIdx < materialLawParams_.size());
00315 return materialLawParams_[elemIdx];
00316 }
00317
00318 template <class FluidState>
00319 void updateHysteresis(const FluidState& fluidState, unsigned elemIdx)
00320 {
00321 if (!enableHysteresis())
00322 return;
00323
00324 auto threePhaseParams = materialLawParams_[elemIdx];
00325 MaterialLaw::updateHysteresis(*threePhaseParams, fluidState);
00326 }
00327
00328 void oilWaterHysteresisParams(Scalar& pcSwMdc,
00329 Scalar& krnSwMdc,
00330 unsigned elemIdx) const
00331 {
00332 if (!enableHysteresis()) {
00333 OPM_THROW(std::runtime_error, "Cannot get hysteresis parameters if hysteresis not enabled.");
00334 }
00335 const auto& params = materialLawParams(elemIdx);
00336 MaterialLaw::oilWaterHysteresisParams(pcSwMdc, krnSwMdc, params);
00337 }
00338
00339 void setOilWaterHysteresisParams(const Scalar& pcSwMdc,
00340 const Scalar& krnSwMdc,
00341 unsigned elemIdx)
00342 {
00343 if (!enableHysteresis()) {
00344 OPM_THROW(std::runtime_error, "Cannot set hysteresis parameters if hysteresis not enabled.");
00345 }
00346 auto& params = materialLawParams(elemIdx);
00347 MaterialLaw::setOilWaterHysteresisParams(pcSwMdc, krnSwMdc, params);
00348 }
00349
00350 void gasOilHysteresisParams(Scalar& pcSwMdc,
00351 Scalar& krnSwMdc,
00352 unsigned elemIdx) const
00353 {
00354 if (!enableHysteresis()) {
00355 OPM_THROW(std::runtime_error, "Cannot get hysteresis parameters if hysteresis not enabled.");
00356 }
00357 const auto& params = materialLawParams(elemIdx);
00358 MaterialLaw::gasOilHysteresisParams(pcSwMdc, krnSwMdc, params);
00359 }
00360
00361 void setGasOilHysteresisParams(const Scalar& pcSwMdc,
00362 const Scalar& krnSwMdc,
00363 unsigned elemIdx)
00364 {
00365 if (!enableHysteresis()) {
00366 OPM_THROW(std::runtime_error, "Cannot set hysteresis parameters if hysteresis not enabled.");
00367 }
00368 auto& params = materialLawParams(elemIdx);
00369 MaterialLaw::setGasOilHysteresisParams(pcSwMdc, krnSwMdc, params);
00370 }
00371
00372 EclEpsScalingPoints<Scalar>& oilWaterScaledEpsPointsDrainage(unsigned elemIdx)
00373 {
00374 auto& materialParams = *materialLawParams_[elemIdx];
00375 switch (materialParams.approach()) {
00376 case EclStone1Approach: {
00377 auto& realParams = materialParams.template getRealParams<Opm::EclStone1Approach>();
00378 return realParams.oilWaterParams().drainageParams().scaledPoints();
00379 }
00380
00381 case EclStone2Approach: {
00382 auto& realParams = materialParams.template getRealParams<Opm::EclStone2Approach>();
00383 return realParams.oilWaterParams().drainageParams().scaledPoints();
00384 }
00385
00386 case EclDefaultApproach: {
00387 auto& realParams = materialParams.template getRealParams<Opm::EclDefaultApproach>();
00388 return realParams.oilWaterParams().drainageParams().scaledPoints();
00389 }
00390
00391 case EclTwoPhaseApproach: {
00392 auto& realParams = materialParams.template getRealParams<Opm::EclTwoPhaseApproach>();
00393 return realParams.oilWaterParams().drainageParams().scaledPoints();
00394 }
00395 default:
00396 OPM_THROW(std::logic_error, "Enum value for material approach unknown!");
00397 }
00398 }
00399
00400 const Opm::EclEpsScalingPointsInfo<Scalar>& oilWaterScaledEpsInfoDrainage(size_t elemIdx) const
00401 {
00402 return *oilWaterScaledEpsInfoDrainage_[elemIdx];
00403 }
00404
00405 std::shared_ptr<EclEpsScalingPointsInfo<Scalar> >& oilWaterScaledEpsInfoDrainagePointerReferenceHack(unsigned elemIdx)
00406 {
00407 return oilWaterScaledEpsInfoDrainage_[elemIdx];
00408 }
00409 private:
00410 void readGlobalEpsOptions_(const Opm::Deck& deck, const Opm::EclipseState& eclState)
00411 {
00412 oilWaterEclEpsConfig_ = std::make_shared<Opm::EclEpsConfig>();
00413 oilWaterEclEpsConfig_-> initFromDeck(deck, eclState, Opm::EclOilWaterSystem);
00414
00415 enableEndPointScaling_ = deck.hasKeyword("ENDSCALE");
00416 }
00417
00418 void readGlobalHysteresisOptions_(const Opm::Deck& deck)
00419 {
00420 hysteresisConfig_ = std::make_shared<Opm::EclHysteresisConfig>();
00421 hysteresisConfig_->initFromDeck(deck);
00422 }
00423
00424 void readGlobalThreePhaseOptions_(const Opm::Deck& deck)
00425 {
00426 bool gasEnabled = deck.hasKeyword("GAS");
00427 bool oilEnabled = deck.hasKeyword("OIL");
00428 bool waterEnabled = deck.hasKeyword("WATER");
00429
00430 int numEnabled =
00431 (gasEnabled?1:0)
00432 + (oilEnabled?1:0)
00433 + (waterEnabled?1:0);
00434
00435 if (numEnabled < 2)
00436 OPM_THROW(std::runtime_error,
00437 "At least two fluid phases must be enabled. (Is: " << numEnabled << ")");
00438
00439 if (numEnabled == 2) {
00440 threePhaseApproach_ = Opm::EclTwoPhaseApproach;
00441 if (!gasEnabled)
00442 twoPhaseApproach_ = Opm::EclTwoPhaseOilWater;
00443 else if (!oilEnabled)
00444 twoPhaseApproach_ = Opm::EclTwoPhaseGasWater;
00445 else if (!waterEnabled)
00446 twoPhaseApproach_ = Opm::EclTwoPhaseGasOil;
00447 }
00448 else {
00449 assert(numEnabled == 3);
00450
00451 threePhaseApproach_ = Opm::EclDefaultApproach;
00452 if (deck.hasKeyword("STONE") || deck.hasKeyword("STONE2"))
00453 threePhaseApproach_ = Opm::EclStone2Approach;
00454 else if (deck.hasKeyword("STONE1"))
00455 threePhaseApproach_ = Opm::EclStone1Approach;
00456 }
00457 }
00458
00459 void initParamsForElements_(const Deck& deck, const EclipseState& eclState,
00460 const std::vector<int>& compressedToCartesianElemIdx,
00461 const std::vector<int>& satnumRegionArray)
00462 {
00463 const size_t numSatRegions = eclState.runspec().tabdims().getNumSatTables();
00464 unsigned numCompressedElems = static_cast<unsigned>(compressedToCartesianElemIdx.size());
00465
00466
00467
00468 auto gasOilConfig = std::make_shared<Opm::EclEpsConfig>();
00469 auto oilWaterConfig = std::make_shared<Opm::EclEpsConfig>();
00470 gasOilConfig->initFromDeck(deck, eclState, Opm::EclGasOilSystem);
00471 oilWaterConfig->initFromDeck(deck, eclState, Opm::EclOilWaterSystem);
00472
00473
00474 gasOilUnscaledPointsVector_.resize(numSatRegions);
00475 oilWaterUnscaledPointsVector_.resize(numSatRegions);
00476 gasOilEffectiveParamVector_.resize(numSatRegions);
00477 oilWaterEffectiveParamVector_.resize(numSatRegions);
00478 for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) {
00479
00480 readGasOilUnscaledPoints_(gasOilUnscaledPointsVector_, gasOilConfig, deck, eclState, satRegionIdx);
00481 readOilWaterUnscaledPoints_(oilWaterUnscaledPointsVector_, oilWaterConfig, deck, eclState, satRegionIdx);
00482
00483
00484 readGasOilEffectiveParameters_(gasOilEffectiveParamVector_, deck, eclState, satRegionIdx);
00485 readOilWaterEffectiveParameters_(oilWaterEffectiveParamVector_, deck, eclState, satRegionIdx);
00486
00487
00488 unscaledEpsInfo_[satRegionIdx].extractUnscaled(deck, eclState, satRegionIdx);
00489
00490 }
00491
00492
00493
00494 GasOilScalingInfoVector gasOilScaledInfoVector(numCompressedElems);
00495 oilWaterScaledEpsInfoDrainage_.resize(numCompressedElems);
00496 GasOilScalingInfoVector gasOilScaledImbInfoVector;
00497 OilWaterScalingInfoVector oilWaterScaledImbInfoVector;
00498
00499 GasOilScalingPointsVector gasOilScaledPointsVector(numCompressedElems);
00500 GasOilScalingPointsVector oilWaterScaledEpsPointsDrainage(numCompressedElems);
00501 GasOilScalingPointsVector gasOilScaledImbPointsVector;
00502 OilWaterScalingPointsVector oilWaterScaledImbPointsVector;
00503
00504 if (enableHysteresis()) {
00505 gasOilScaledImbInfoVector.resize(numCompressedElems);
00506 gasOilScaledImbPointsVector.resize(numCompressedElems);
00507 oilWaterScaledImbInfoVector.resize(numCompressedElems);
00508 oilWaterScaledImbPointsVector.resize(numCompressedElems);
00509 }
00510
00511 EclEpsGridProperties epsGridProperties, epsImbGridProperties;
00512 epsGridProperties.initFromDeck(deck, eclState, false);
00513 if (enableHysteresis())
00514 epsImbGridProperties.initFromDeck(deck, eclState, true);
00515 for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
00516 unsigned cartElemIdx = static_cast<unsigned>(compressedToCartesianElemIdx[elemIdx]);
00517 readGasOilScaledPoints_(gasOilScaledInfoVector,
00518 gasOilScaledPointsVector,
00519 gasOilConfig,
00520 eclState,
00521 epsGridProperties,
00522 elemIdx,
00523 cartElemIdx);
00524 readOilWaterScaledPoints_(oilWaterScaledEpsInfoDrainage_,
00525 oilWaterScaledEpsPointsDrainage,
00526 oilWaterConfig,
00527 eclState,
00528 epsGridProperties,
00529 elemIdx,
00530 cartElemIdx);
00531
00532 if (enableHysteresis()) {
00533 readGasOilScaledPoints_(gasOilScaledImbInfoVector,
00534 gasOilScaledImbPointsVector,
00535 gasOilConfig,
00536 eclState,
00537 epsImbGridProperties,
00538 elemIdx,
00539 cartElemIdx);
00540 readOilWaterScaledPoints_(oilWaterScaledImbInfoVector,
00541 oilWaterScaledImbPointsVector,
00542 oilWaterConfig,
00543 eclState,
00544 epsImbGridProperties,
00545 elemIdx,
00546 cartElemIdx);
00547 }
00548 }
00549
00550
00551 GasOilParamVector gasOilParams(numCompressedElems);
00552 OilWaterParamVector oilWaterParams(numCompressedElems);
00553 GasOilParamVector gasOilImbParams;
00554 OilWaterParamVector oilWaterImbParams;
00555
00556 if (enableHysteresis()) {
00557 gasOilImbParams.resize(numCompressedElems);
00558 oilWaterImbParams.resize(numCompressedElems);
00559 }
00560
00561 bool hasGas = deck.hasKeyword("GAS");
00562 bool hasOil = deck.hasKeyword("OIL");
00563 bool hasWater = deck.hasKeyword("WATER");
00564
00565 const auto& imbnumData = eclState.get3DProperties().getIntGridProperty("IMBNUM").getData();
00566 assert(numCompressedElems == satnumRegionArray.size());
00567 for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
00568 unsigned satRegionIdx = static_cast<unsigned>(satnumRegionArray[elemIdx]);
00569
00570 gasOilParams[elemIdx] = std::make_shared<GasOilTwoPhaseHystParams>();
00571 oilWaterParams[elemIdx] = std::make_shared<OilWaterTwoPhaseHystParams>();
00572
00573 gasOilParams[elemIdx]->setConfig(hysteresisConfig_);
00574 oilWaterParams[elemIdx]->setConfig(hysteresisConfig_);
00575
00576 if (hasGas && hasOil) {
00577 auto gasOilDrainParams = std::make_shared<GasOilEpsTwoPhaseParams>();
00578 gasOilDrainParams->setConfig(gasOilConfig);
00579 gasOilDrainParams->setUnscaledPoints(gasOilUnscaledPointsVector_[satRegionIdx]);
00580 gasOilDrainParams->setScaledPoints(gasOilScaledPointsVector[elemIdx]);
00581 gasOilDrainParams->setEffectiveLawParams(gasOilEffectiveParamVector_[satRegionIdx]);
00582 gasOilDrainParams->finalize();
00583
00584 gasOilParams[elemIdx]->setDrainageParams(gasOilDrainParams,
00585 *gasOilScaledInfoVector[elemIdx],
00586 EclGasOilSystem);
00587 }
00588
00589 if (hasOil && hasWater) {
00590 auto oilWaterDrainParams = std::make_shared<OilWaterEpsTwoPhaseParams>();
00591 oilWaterDrainParams->setConfig(oilWaterConfig);
00592 oilWaterDrainParams->setUnscaledPoints(oilWaterUnscaledPointsVector_[satRegionIdx]);
00593 oilWaterDrainParams->setScaledPoints(oilWaterScaledEpsPointsDrainage[elemIdx]);
00594 oilWaterDrainParams->setEffectiveLawParams(oilWaterEffectiveParamVector_[satRegionIdx]);
00595 oilWaterDrainParams->finalize();
00596
00597 oilWaterParams[elemIdx]->setDrainageParams(oilWaterDrainParams,
00598 *oilWaterScaledEpsInfoDrainage_[elemIdx],
00599 EclOilWaterSystem);
00600 }
00601
00602 if (enableHysteresis()) {
00603 unsigned imbRegionIdx = static_cast<unsigned>(imbnumData[elemIdx]) - 1;
00604
00605 if (hasGas && hasOil) {
00606 auto gasOilImbParamsHyst = std::make_shared<GasOilEpsTwoPhaseParams>();
00607 gasOilImbParamsHyst->setConfig(gasOilConfig);
00608 gasOilImbParamsHyst->setUnscaledPoints(gasOilUnscaledPointsVector_[imbRegionIdx]);
00609 gasOilImbParamsHyst->setScaledPoints(gasOilScaledImbPointsVector[elemIdx]);
00610 gasOilImbParamsHyst->setEffectiveLawParams(gasOilEffectiveParamVector_[imbRegionIdx]);
00611 gasOilImbParamsHyst->finalize();
00612
00613 gasOilParams[elemIdx]->setImbibitionParams(gasOilImbParamsHyst,
00614 *gasOilScaledImbInfoVector[elemIdx],
00615 EclGasOilSystem);
00616 }
00617
00618 if (hasOil && hasWater) {
00619 auto oilWaterImbParamsHyst = std::make_shared<OilWaterEpsTwoPhaseParams>();
00620 oilWaterImbParamsHyst->setConfig(oilWaterConfig);
00621 oilWaterImbParamsHyst->setUnscaledPoints(oilWaterUnscaledPointsVector_[imbRegionIdx]);
00622 oilWaterImbParamsHyst->setScaledPoints(oilWaterScaledImbPointsVector[elemIdx]);
00623 oilWaterImbParamsHyst->setEffectiveLawParams(oilWaterEffectiveParamVector_[imbRegionIdx]);
00624 oilWaterImbParamsHyst->finalize();
00625
00626 oilWaterParams[elemIdx]->setImbibitionParams(oilWaterImbParamsHyst,
00627 *gasOilScaledImbInfoVector[elemIdx],
00628 EclGasOilSystem);
00629 }
00630 }
00631
00632 if (hasGas && hasOil)
00633 gasOilParams[elemIdx]->finalize();
00634
00635 if (hasOil && hasWater)
00636 oilWaterParams[elemIdx]->finalize();
00637 }
00638
00639
00640 materialLawParams_.resize(numCompressedElems);
00641 for (unsigned elemIdx = 0; elemIdx < numCompressedElems; ++elemIdx) {
00642 materialLawParams_[elemIdx] = std::make_shared<MaterialLawParams>();
00643 unsigned satRegionIdx = static_cast<unsigned>(satnumRegionArray[elemIdx]);
00644
00645 initThreePhaseParams_(deck,
00646 eclState,
00647 *materialLawParams_[elemIdx],
00648 satRegionIdx,
00649 *oilWaterScaledEpsInfoDrainage_[elemIdx],
00650 oilWaterParams[elemIdx],
00651 gasOilParams[elemIdx]);
00652
00653 materialLawParams_[elemIdx]->finalize();
00654 }
00655 }
00656
00657
00658
00659
00660
00661 enum SaturationFunctionFamily {
00662 noFamily,
00663 FamilyI,
00664 FamilyII
00665 };
00666
00667 SaturationFunctionFamily getSaturationFunctionFamily(const Opm::Deck& deck, const Opm::EclipseState& eclState) const
00668 {
00669 const auto& tableManager = eclState.getTableManager();
00670 const TableContainer& swofTables = tableManager.getSwofTables();
00671 const TableContainer& slgofTables= tableManager.getSlgofTables();
00672 const TableContainer& sgofTables = tableManager.getSgofTables();
00673 const TableContainer& swfnTables = tableManager.getSwfnTables();
00674 const TableContainer& sgfnTables = tableManager.getSgfnTables();
00675 const TableContainer& sof3Tables = tableManager.getSof3Tables();
00676
00677 bool hasGas = deck.hasKeyword("GAS");
00678 bool hasOil = deck.hasKeyword("OIL");
00679 bool hasWater = deck.hasKeyword("WATER");
00680
00681 bool family1 = false;
00682 bool family2 = false;
00683 if (!hasGas) {
00684
00685 family1 = !swofTables.empty();
00686 if (!family1)
00687 throw std::runtime_error("only SWOF is supporeted for oil-water two-phase simulations");
00688
00689 }
00690 else if (!hasWater) {
00691
00692 family1 = !sgofTables.empty();
00693 if (!family1)
00694 throw std::runtime_error("only SGOF is supporeted for oil-gas two-phase simulations");
00695
00696 }
00697 else if (!hasOil) {
00698
00699 throw std::runtime_error("water-gas two-phase simulations are currently not supported");
00700 }
00701 else {
00702
00703 family1 = (!sgofTables.empty() || !slgofTables.empty()) && !swofTables.empty();
00704 family2 = !swfnTables.empty() && !sgfnTables.empty() && !sof3Tables.empty();
00705 }
00706
00707 if (family1 && family2) {
00708 throw std::invalid_argument("Saturation families should not be mixed \n"
00709 "Use either SGOF and SWOF or SGFN, SWFN and SOF3");
00710 }
00711
00712 if (!family1 && !family2) {
00713 throw std::invalid_argument("Saturations function must be specified using either "
00714 "family 1 or family 2 keywords \n"
00715 "Use either SGOF and SWOF or SGFN, SWFN and SOF3" );
00716 }
00717
00718 if (family1 && !family2)
00719 return SaturationFunctionFamily::FamilyI;
00720 else if (family2 && !family1)
00721 return SaturationFunctionFamily::FamilyII;
00722 return SaturationFunctionFamily::noFamily;
00723 }
00724
00725 template <class Container>
00726 void readGasOilEffectiveParameters_(Container& dest,
00727 const Opm::Deck& deck,
00728 const Opm::EclipseState& eclState,
00729 unsigned satRegionIdx)
00730 {
00731 bool hasGas = deck.hasKeyword("GAS");
00732 bool hasOil = deck.hasKeyword("OIL");
00733
00734 if (!hasGas || !hasOil)
00735
00736 return;
00737
00738 dest[satRegionIdx] = std::make_shared<GasOilEffectiveTwoPhaseParams>();
00739
00740 auto& effParams = *dest[satRegionIdx];
00741
00742
00743
00744 Scalar Swco = unscaledEpsInfo_[satRegionIdx].Swl;
00745 const auto& tableManager = eclState.getTableManager();
00746
00747 switch (getSaturationFunctionFamily(deck, eclState)) {
00748 case FamilyI:
00749 {
00750 const TableContainer& sgofTables = tableManager.getSgofTables();
00751 const TableContainer& slgofTables = tableManager.getSlgofTables();
00752 if (!sgofTables.empty())
00753 readGasOilEffectiveParametersSgof_(effParams,
00754 Swco,
00755 sgofTables.getTable<SgofTable>(satRegionIdx));
00756 else if (!slgofTables.empty())
00757 readGasOilEffectiveParametersSlgof_(effParams,
00758 Swco,
00759 slgofTables.getTable<SlgofTable>(satRegionIdx));
00760 break;
00761 }
00762
00763 case FamilyII:
00764 {
00765 const Sof3Table& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>( satRegionIdx );
00766 const SgfnTable& sgfnTable = tableManager.getSgfnTables().getTable<SgfnTable>( satRegionIdx );
00767 readGasOilEffectiveParametersFamily2_(effParams,
00768 Swco,
00769 sof3Table,
00770 sgfnTable);
00771 break;
00772 }
00773
00774
00775 case noFamily:
00776 throw std::domain_error("No valid saturation keyword family specified");
00777
00778 }
00779 }
00780
00781 void readGasOilEffectiveParametersSgof_(GasOilEffectiveTwoPhaseParams& effParams,
00782 Scalar Swco,
00783 const Opm::SgofTable& sgofTable)
00784 {
00785
00786 std::vector<double> SoSamples(sgofTable.numRows());
00787 std::vector<double> SoKroSamples(sgofTable.numRows());
00788 for (size_t sampleIdx = 0; sampleIdx < sgofTable.numRows(); ++ sampleIdx) {
00789 SoSamples[sampleIdx] = 1 - sgofTable.get("SG", sampleIdx);
00790 SoKroSamples[sampleIdx] = SoSamples[sampleIdx] - Swco;
00791 }
00792
00793 effParams.setKrwSamples(SoKroSamples, sgofTable.getColumn("KROG").vectorCopy());
00794 effParams.setKrnSamples(SoSamples, sgofTable.getColumn("KRG").vectorCopy());
00795 effParams.setPcnwSamples(SoSamples, sgofTable.getColumn("PCOG").vectorCopy());
00796 effParams.finalize();
00797 }
00798
00799 void readGasOilEffectiveParametersSlgof_(GasOilEffectiveTwoPhaseParams& effParams,
00800 Scalar Swco,
00801 const Opm::SlgofTable& slgofTable)
00802 {
00803
00804 std::vector<double> SoSamples(slgofTable.numRows());
00805 std::vector<double> SoKroSamples(slgofTable.numRows());
00806 for (size_t sampleIdx = 0; sampleIdx < slgofTable.numRows(); ++ sampleIdx) {
00807 SoSamples[sampleIdx] = slgofTable.get("SL", sampleIdx);
00808 SoKroSamples[sampleIdx] = slgofTable.get("SL", sampleIdx) - Swco;
00809 }
00810
00811 effParams.setKrwSamples(SoKroSamples, slgofTable.getColumn("KROG").vectorCopy());
00812 effParams.setKrnSamples(SoSamples, slgofTable.getColumn("KRG").vectorCopy());
00813 effParams.setPcnwSamples(SoSamples, slgofTable.getColumn("PCOG").vectorCopy());
00814 effParams.finalize();
00815 }
00816
00817 void readGasOilEffectiveParametersFamily2_(GasOilEffectiveTwoPhaseParams& effParams,
00818 Scalar ,
00819 const Opm::Sof3Table& sof3Table,
00820 const Opm::SgfnTable& sgfnTable)
00821 {
00822
00823 std::vector<double> SoSamples(sgfnTable.numRows());
00824 std::vector<double> SoColumn = sof3Table.getColumn("SO").vectorCopy();
00825 for (size_t sampleIdx = 0; sampleIdx < sgfnTable.numRows(); ++ sampleIdx) {
00826 SoSamples[sampleIdx] = 1 - sgfnTable.get("SG", sampleIdx);
00827 }
00828
00829 effParams.setKrwSamples(SoColumn, sof3Table.getColumn("KROG").vectorCopy());
00830 effParams.setKrnSamples(SoSamples, sgfnTable.getColumn("KRG").vectorCopy());
00831 effParams.setPcnwSamples(SoSamples, sgfnTable.getColumn("PCOG").vectorCopy());
00832 effParams.finalize();
00833 }
00834
00835 template <class Container>
00836 void readOilWaterEffectiveParameters_(Container& dest,
00837 const Opm::Deck& deck,
00838 const Opm::EclipseState& eclState,
00839 unsigned satRegionIdx)
00840 {
00841 bool hasWater = deck.hasKeyword("WATER");
00842 bool hasOil = deck.hasKeyword("OIL");
00843
00844 if (!hasOil || !hasWater)
00845
00846 return;
00847
00848 dest[satRegionIdx] = std::make_shared<OilWaterEffectiveTwoPhaseParams>();
00849
00850 const auto& tableManager = eclState.getTableManager();
00851 auto& effParams = *dest[satRegionIdx];
00852
00853 switch (getSaturationFunctionFamily(deck, eclState)) {
00854 case FamilyI: {
00855 const auto& swofTable = tableManager.getSwofTables().getTable<SwofTable>(satRegionIdx);
00856 std::vector<double> SwColumn = swofTable.getColumn("SW").vectorCopy();
00857
00858 effParams.setKrwSamples(SwColumn, swofTable.getColumn("KRW").vectorCopy());
00859 effParams.setKrnSamples(SwColumn, swofTable.getColumn("KROW").vectorCopy());
00860 effParams.setPcnwSamples(SwColumn, swofTable.getColumn("PCOW").vectorCopy());
00861 effParams.finalize();
00862 break;
00863 }
00864 case FamilyII:
00865 {
00866 const auto& swfnTable = tableManager.getSwfnTables().getTable<SwfnTable>(satRegionIdx);
00867 const auto& sof3Table = tableManager.getSof3Tables().getTable<Sof3Table>(satRegionIdx);
00868 std::vector<double> SwColumn = swfnTable.getColumn("SW").vectorCopy();
00869
00870
00871 std::vector<double> SwSamples(sof3Table.numRows());
00872 for (size_t sampleIdx = 0; sampleIdx < sof3Table.numRows(); ++ sampleIdx)
00873 SwSamples[sampleIdx] = 1 - sof3Table.get("SO", sampleIdx);
00874
00875 effParams.setKrwSamples(SwColumn, swfnTable.getColumn("KRW").vectorCopy());
00876 effParams.setKrnSamples(SwSamples, sof3Table.getColumn("KROW").vectorCopy());
00877 effParams.setPcnwSamples(SwColumn, swfnTable.getColumn("PCOW").vectorCopy());
00878 effParams.finalize();
00879 break;
00880 }
00881
00882 case noFamily:
00883
00884 throw std::domain_error("No valid saturation keyword family specified");
00885
00886 }
00887 }
00888
00889
00890 template <class Container>
00891 void readGasOilUnscaledPoints_(Container& dest,
00892 std::shared_ptr<EclEpsConfig> config,
00893 const Opm::Deck& deck,
00894 const Opm::EclipseState& ,
00895 unsigned satRegionIdx)
00896 {
00897 bool hasGas = deck.hasKeyword("GAS");
00898 bool hasOil = deck.hasKeyword("OIL");
00899
00900 if (!hasGas || !hasOil)
00901
00902 return;
00903
00904 dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
00905 dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], *config, EclGasOilSystem);
00906 }
00907
00908 template <class Container>
00909 void readOilWaterUnscaledPoints_(Container& dest,
00910 std::shared_ptr<EclEpsConfig> config,
00911 const Opm::Deck& deck,
00912 const Opm::EclipseState& ,
00913 unsigned satRegionIdx)
00914 {
00915 bool hasWater = deck.hasKeyword("WATER");
00916 bool hasOil = deck.hasKeyword("OIL");
00917
00918 if (!hasOil || !hasWater)
00919
00920 return;
00921
00922 dest[satRegionIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
00923 dest[satRegionIdx]->init(unscaledEpsInfo_[satRegionIdx], *config, EclOilWaterSystem);
00924 }
00925
00926 template <class InfoContainer, class PointsContainer>
00927 void readGasOilScaledPoints_(InfoContainer& destInfo,
00928 PointsContainer& destPoints,
00929 std::shared_ptr<EclEpsConfig> config,
00930 const Opm::EclipseState& eclState,
00931 const EclEpsGridProperties& epsGridProperties,
00932 unsigned elemIdx,
00933 unsigned cartElemIdx)
00934 {
00935 unsigned satRegionIdx = static_cast<unsigned>((*epsGridProperties.satnum)[cartElemIdx]) - 1;
00936
00937 destInfo[elemIdx] = std::make_shared<EclEpsScalingPointsInfo<Scalar> >(unscaledEpsInfo_[satRegionIdx]);
00938 destInfo[elemIdx]->extractScaled(eclState, epsGridProperties, cartElemIdx);
00939
00940 destPoints[elemIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
00941 destPoints[elemIdx]->init(*destInfo[elemIdx], *config, EclGasOilSystem);
00942 }
00943
00944 template <class InfoContainer, class PointsContainer>
00945 void readOilWaterScaledPoints_(InfoContainer& destInfo,
00946 PointsContainer& destPoints,
00947 std::shared_ptr<EclEpsConfig> config,
00948 const Opm::EclipseState& eclState,
00949 const EclEpsGridProperties& epsGridProperties,
00950 unsigned elemIdx,
00951 unsigned cartElemIdx)
00952 {
00953 unsigned satRegionIdx = static_cast<unsigned>((*epsGridProperties.satnum)[cartElemIdx]) - 1;
00954
00955 destInfo[elemIdx] = std::make_shared<EclEpsScalingPointsInfo<Scalar> >(unscaledEpsInfo_[satRegionIdx]);
00956 destInfo[elemIdx]->extractScaled(eclState, epsGridProperties, cartElemIdx);
00957
00958 destPoints[elemIdx] = std::make_shared<EclEpsScalingPoints<Scalar> >();
00959 destPoints[elemIdx]->init(*destInfo[elemIdx], *config, EclOilWaterSystem);
00960 }
00961
00962 void initThreePhaseParams_(const Opm::Deck& deck,
00963 const Opm::EclipseState& ,
00964 MaterialLawParams& materialParams,
00965 unsigned satRegionIdx,
00966 const EclEpsScalingPointsInfo<Scalar>& epsInfo,
00967 std::shared_ptr<OilWaterTwoPhaseHystParams> oilWaterParams,
00968 std::shared_ptr<GasOilTwoPhaseHystParams> gasOilParams)
00969 {
00970 materialParams.setApproach(threePhaseApproach_);
00971
00972 switch (materialParams.approach()) {
00973 case EclStone1Approach: {
00974 auto& realParams = materialParams.template getRealParams<Opm::EclStone1Approach>();
00975 realParams.setGasOilParams(gasOilParams);
00976 realParams.setOilWaterParams(oilWaterParams);
00977 realParams.setSwl(epsInfo.Swl);
00978
00979 if (deck.hasKeyword("STONE1EX")) {
00980 Scalar eta =
00981 deck.getKeyword("STONE1EX").getRecord(satRegionIdx).getItem(0).getSIDouble(0);
00982 realParams.setEta(eta);
00983 }
00984 else
00985 realParams.setEta(1.0);
00986 realParams.finalize();
00987 break;
00988 }
00989
00990 case EclStone2Approach: {
00991 auto& realParams = materialParams.template getRealParams<Opm::EclStone2Approach>();
00992 realParams.setGasOilParams(gasOilParams);
00993 realParams.setOilWaterParams(oilWaterParams);
00994 realParams.setSwl(epsInfo.Swl);
00995 realParams.finalize();
00996 break;
00997 }
00998
00999 case EclDefaultApproach: {
01000 auto& realParams = materialParams.template getRealParams<Opm::EclDefaultApproach>();
01001 realParams.setGasOilParams(gasOilParams);
01002 realParams.setOilWaterParams(oilWaterParams);
01003 realParams.setSwl(epsInfo.Swl);
01004 realParams.finalize();
01005 break;
01006 }
01007
01008 case EclTwoPhaseApproach: {
01009 auto& realParams = materialParams.template getRealParams<Opm::EclTwoPhaseApproach>();
01010 realParams.setGasOilParams(gasOilParams);
01011 realParams.setOilWaterParams(oilWaterParams);
01012 realParams.setApproach(twoPhaseApproach_);
01013 realParams.finalize();
01014 break;
01015 }
01016 }
01017 }
01018
01019 bool enableEndPointScaling_;
01020 std::shared_ptr<EclHysteresisConfig> hysteresisConfig_;
01021
01022 std::shared_ptr<EclEpsConfig> oilWaterEclEpsConfig_;
01023 std::vector<Opm::EclEpsScalingPointsInfo<Scalar>> unscaledEpsInfo_;
01024 OilWaterScalingInfoVector oilWaterScaledEpsInfoDrainage_;
01025
01026 GasOilScalingPointsVector gasOilUnscaledPointsVector_;
01027 OilWaterScalingPointsVector oilWaterUnscaledPointsVector_;
01028 GasOilEffectiveParamVector gasOilEffectiveParamVector_;
01029 OilWaterEffectiveParamVector oilWaterEffectiveParamVector_;
01030
01031 Opm::EclMultiplexerApproach threePhaseApproach_;
01032
01033 enum EclTwoPhaseApproach twoPhaseApproach_;
01034
01035 std::vector<std::shared_ptr<MaterialLawParams> > materialLawParams_;
01036
01037 std::vector<int> satnumRegionArray_;
01038 };
01039 }
01040
01041 #endif