From 4f84ca4716081e45c3a78fb67d8b13660eb1a774 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 25 Sep 2020 15:10:39 +0200 Subject: [PATCH] Add support for ROCKTAB --- ebos/eclproblem.hh | 76 ++++++++++++++++++++++++++++++++-------------- 1 file changed, 54 insertions(+), 22 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index dfcb1c63f..fadb168cc 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -81,6 +81,7 @@ #include #include #include +#include #include #include @@ -92,6 +93,7 @@ #include #include #include +#include #include #include @@ -656,6 +658,7 @@ class EclProblem : public GetPropType typedef typename GridView::template Codim<0>::Iterator ElementIterator; typedef Opm::UniformXTabulated2DFunction TabulatedTwoDFunction; + typedef Opm::Tabulated1DFunction TabulatedFunction; struct RockParams { Scalar referencePressure; @@ -1121,7 +1124,7 @@ public: * using a pore volume multiplier (i.e., poreVolumeMultiplier() != 1.0) */ bool recycleFirstIterationStorage() const - { return !drsdtActive_() && !drvdtActive_() && rockCompPoroMult_.empty(); } + { return !drsdtActive_() && !drvdtActive_() && rockCompPoroMultWc_.empty() && rockCompPoroMult_.empty(); } /*! * \brief Called by the simulator before each Newton-Raphson iteration. @@ -2035,7 +2038,8 @@ public: template LhsEval rockCompPoroMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const { - if (rockCompPoroMult_.empty()) + + if (rockCompPoroMult_.empty() && rockCompPoroMultWc_.empty()) return 1.0; unsigned tableIdx = 0; @@ -2043,10 +2047,7 @@ public: tableIdx = rockTableIdx_[elementIdx]; const auto& fs = intQuants.fluidState(); - LhsEval SwMax = Opm::max(Opm::decay(fs.saturation(waterPhaseIdx)), maxWaterSaturation_[elementIdx]); - LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx); LhsEval effectiveOilPressure = Opm::decay(fs.pressure(oilPhaseIdx)); - if (!minOilPressure_.empty()) // The pore space change is irreversible effectiveOilPressure = @@ -2056,7 +2057,17 @@ public: if (!overburdenPressure_.empty()) effectiveOilPressure -= overburdenPressure_[elementIdx]; - return rockCompPoroMult_[tableIdx].eval(effectiveOilPressure, SwDeltaMax, /*extrapolation=*/true); + + if (!rockCompPoroMult_.empty()) { + return rockCompPoroMult_[tableIdx].eval(effectiveOilPressure, /*extrapolation=*/true); + } + + // water compaction + assert(!rockCompPoroMultWc_.empty()); + LhsEval SwMax = Opm::max(Opm::decay(fs.saturation(waterPhaseIdx)), maxWaterSaturation_[elementIdx]); + LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx); + + return rockCompPoroMultWc_[tableIdx].eval(effectiveOilPressure, SwDeltaMax, /*extrapolation=*/true); } /*! @@ -2067,7 +2078,7 @@ public: template LhsEval rockCompTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const { - if (rockCompTransMult_.empty()) + if (rockCompTransMult_.empty() && rockCompTransMultWc_.empty()) return 1.0; unsigned tableIdx = 0; @@ -2075,8 +2086,6 @@ public: tableIdx = rockTableIdx_[elementIdx]; const auto& fs = intQuants.fluidState(); - LhsEval SwMax = Opm::max(Opm::decay(fs.saturation(waterPhaseIdx)), maxWaterSaturation_[elementIdx]); - LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx); LhsEval effectiveOilPressure = Opm::decay(fs.pressure(oilPhaseIdx)); if (!minOilPressure_.empty()) @@ -2088,7 +2097,15 @@ public: if (overburdenPressure_.size() > 0) effectiveOilPressure -= overburdenPressure_[elementIdx]; - return rockCompTransMult_[tableIdx].eval(effectiveOilPressure, SwDeltaMax, /*extrapolation=*/true); + if (!rockCompTransMult_.empty()) + return rockCompTransMult_[tableIdx].eval(effectiveOilPressure, /*extrapolation=*/true); + + // water compaction + assert(!rockCompTransMultWc_.empty()); + LhsEval SwMax = Opm::max(Opm::decay(fs.saturation(waterPhaseIdx)), maxWaterSaturation_[elementIdx]); + LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx); + + return rockCompTransMultWc_[tableIdx].eval(effectiveOilPressure, SwDeltaMax, /*extrapolation=*/true); } /*! @@ -2392,8 +2409,8 @@ private: readRockCompactionParameters_(); unsigned numElem = vanguard.gridView().size(0); - rockTableIdx_.resize(numElem); if (eclState.fieldProps().has_int(rock_config.rocknum_property())) { + rockTableIdx_.resize(numElem); const auto& num = eclState.fieldProps().get_int(rock_config.rocknum_property()); for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) { rockTableIdx_[elemIdx] = num[elemIdx] - 1; @@ -2451,10 +2468,23 @@ private: size_t numRocktabTables = rock_config.num_rock_tables(); bool waterCompaction = rock_config.water_compaction(); - if (!waterCompaction) - throw std::runtime_error("Only water compatction allowed"); + if (!waterCompaction) { + const auto& rocktabTables = eclState.getTableManager().getRocktabTables(); + if (rocktabTables.size() != numRocktabTables) + throw std::runtime_error("ROCKCOMP is activated." + std::to_string(numRocktabTables) + +" ROCKTAB tables is expected, but " + std::to_string(rocktabTables.size()) +" is provided"); - if (waterCompaction) { + rockCompPoroMult_.resize(numRocktabTables); + rockCompTransMult_.resize(numRocktabTables); + for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { + const auto& rocktabTable = rocktabTables.template getTable(regionIdx); + const auto& pressureColumn = rocktabTable.getPressureColumn(); + const auto& poroColumn = rocktabTable.getPoreVolumeMultiplierColumn(); + const auto& transColumn = rocktabTable.getTransmissibilityMultiplierColumn(); + rockCompPoroMult_[regionIdx].setXYContainers(pressureColumn, poroColumn); + rockCompTransMult_[regionIdx].setXYContainers(pressureColumn, transColumn); + } + } else { const auto& rock2dTables = eclState.getTableManager().getRock2dTables(); const auto& rock2dtrTables = eclState.getTableManager().getRock2dtrTables(); const auto& rockwnodTables = eclState.getTableManager().getRockwnodTables(); @@ -2468,7 +2498,7 @@ private: throw std::runtime_error("Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables) +" ROCKWNOD tables is expected, but " + std::to_string(rockwnodTables.size()) +" is provided"); //TODO check size match - rockCompPoroMult_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical)); + rockCompPoroMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical)); for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { const Opm::RockwnodTable& rockwnodTable = rockwnodTables.template getTable(regionIdx); const auto& rock2dTable = rock2dTables[regionIdx]; @@ -2477,16 +2507,16 @@ private: throw std::runtime_error("Number of entries in ROCKWNOD and ROCK2D needs to match."); for (size_t xIdx = 0; xIdx < rock2dTable.size(); ++xIdx) { - rockCompPoroMult_[regionIdx].appendXPos(rock2dTable.getPressureValue(xIdx)); + rockCompPoroMultWc_[regionIdx].appendXPos(rock2dTable.getPressureValue(xIdx)); for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx) - rockCompPoroMult_[regionIdx].appendSamplePoint(xIdx, + rockCompPoroMultWc_[regionIdx].appendSamplePoint(xIdx, rockwnodTable.getSaturationColumn()[yIdx], rock2dTable.getPvmultValue(xIdx, yIdx)); } } if (!rock2dtrTables.empty()) { - rockCompTransMult_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical)); + rockCompTransMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical)); for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { const Opm::RockwnodTable& rockwnodTable = rockwnodTables.template getTable(regionIdx); const auto& rock2dtrTable = rock2dtrTables[regionIdx]; @@ -2495,9 +2525,9 @@ private: throw std::runtime_error("Number of entries in ROCKWNOD and ROCK2DTR needs to match."); for (size_t xIdx = 0; xIdx < rock2dtrTable.size(); ++xIdx) { - rockCompTransMult_[regionIdx].appendXPos(rock2dtrTable.getPressureValue(xIdx)); + rockCompTransMultWc_[regionIdx].appendXPos(rock2dtrTable.getPressureValue(xIdx)); for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx) - rockCompTransMult_[regionIdx].appendSamplePoint(xIdx, + rockCompTransMultWc_[regionIdx].appendSamplePoint(xIdx, rockwnodTable.getSaturationColumn()[yIdx], rock2dtrTable.getTransMultValue(xIdx, yIdx)); } @@ -3280,8 +3310,10 @@ private: std::vector overburdenPressure_; std::vector minOilPressure_; - std::vector rockCompPoroMult_; - std::vector rockCompTransMult_; + std::vector rockCompPoroMultWc_; + std::vector rockCompTransMultWc_; + std::vector rockCompPoroMult_; + std::vector rockCompTransMult_; bool enableDriftCompensation_; GlobalEqVector drift_;