diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index cb5368a0a..399efa8e4 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -119,6 +119,7 @@ class EclTransExtensiveQuantities enum { numPhases = FluidSystem::numPhases }; enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + enum { enableExperiments = GET_PROP_VALUE(TypeTag, EnableExperiments) }; typedef Opm::MathToolbox Toolbox; typedef Dune::FieldVector DimVector; @@ -343,12 +344,18 @@ protected: // does not matter, though. unsigned upstreamIdx = upstreamIndex_(phaseIdx); const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx); + + // TODO: should the rock compaction transmissibility multiplier be upstreamed + // or averaged? all fluids should see the same compaction?! + const Evaluation& transMult = + problem.template rockCompTransMultiplier(up, stencil.globalSpaceIndex(upstreamIdx)); + if (upstreamIdx == interiorDofIdx_) volumeFlux_[phaseIdx] = - pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-trans/faceArea); + pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*transMult*(-trans/faceArea); else volumeFlux_[phaseIdx] = - pressureDifference_[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*(-trans/faceArea)); + pressureDifference_[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*Toolbox::value(transMult)*(-trans/faceArea)); } } @@ -427,14 +434,20 @@ protected: // only works for the element centered finite volume method. for ebos this // does not matter, though. unsigned upstreamIdx = upstreamIndex_(phaseIdx); - if (upstreamIdx == interiorDofIdx_) { - const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx); - volumeFlux_[phaseIdx] = - pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-trans/faceArea); + const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx); - if (enableSolvent && phaseIdx == gasPhaseIdx) { - asImp_().setSolventVolumeFlux( pressureDifference_[phaseIdx]*up.solventMobility()*(-trans/faceArea)); - } + Evaluation transModified = trans; + if (enableExperiments) { + // deal with water induced rock compaction + transModified *= problem.template rockCompTransMultiplier(up, stencil.globalSpaceIndex(upstreamIdx)); + } + + if (upstreamIdx == interiorDofIdx_) { + volumeFlux_[phaseIdx] = + pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea); + + if (enableSolvent && phaseIdx == gasPhaseIdx) + asImp_().setSolventVolumeFlux( pressureDifference_[phaseIdx]*up.solventMobility()*(-transModified/faceArea)); } else { // compute the phase mobility using the material law parameters of the @@ -448,12 +461,11 @@ protected: const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx); volumeFlux_[phaseIdx] = - pressureDifference_[phaseIdx]*mob*(-trans/faceArea); + pressureDifference_[phaseIdx]*mob*(-transModified/faceArea); // Solvent inflow is not yet supported - if (enableSolvent && phaseIdx == gasPhaseIdx) { - asImp_().setSolventVolumeFlux( 0.0 ); - } + if (enableSolvent && phaseIdx == gasPhaseIdx) + asImp_().setSolventVolumeFlux(0.0); } } } diff --git a/ebos/eclnewtonmethod.hh b/ebos/eclnewtonmethod.hh index 5471387ed..ba49abdb2 100644 --- a/ebos/eclnewtonmethod.hh +++ b/ebos/eclnewtonmethod.hh @@ -160,7 +160,7 @@ public: const auto& r = currentResidual[dofIdx]; Scalar pvValue = - this->simulator_.problem().porosity(dofIdx) + this->simulator_.problem().referencePorosity(dofIdx, /*timeIdx=*/0) * this->model().dofTotalVolume(dofIdx); sumPv += pvValue; bool cnvViolated = false; diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 3ecebbe94..5dc60488e 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -335,6 +335,16 @@ public: } } + // ROCKC + if (rstKeywords["ROCKC"] > 0) { + rstKeywords["ROCKC"] = 0; + rockCompPorvMultiplier_.resize(bufferSize, 0.0); + rockCompTransMultiplier_.resize(bufferSize, 0.0); + swMax_.resize(bufferSize, 0.0); + minimumOilPressure_.resize(bufferSize, 0.0); + overburdenPressure_.resize(bufferSize, 0.0); + } + //Warn for any unhandled keyword if (log) { for (auto& keyValue: rstKeywords) { @@ -367,6 +377,7 @@ public: if (!std::is_same >::value) return; + const auto& problem = elemCtx.simulator().problem(); for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); const auto& fs = intQuants.fluidState(); @@ -495,9 +506,30 @@ public: } if (soMax_.size() > 0) - soMax_[globalDofIdx] = elemCtx.simulator().problem().maxOilSaturation(globalDofIdx); + soMax_[globalDofIdx] = + std::max(Opm::getValue(fs.saturation(oilPhaseIdx)), + problem.maxOilSaturation(globalDofIdx)); - const auto& matLawManager = elemCtx.simulator().problem().materialLawManager(); + if (swMax_.size() > 0) + swMax_[globalDofIdx] = + std::max(Opm::getValue(fs.saturation(waterPhaseIdx)), + problem.maxWaterSaturation(globalDofIdx)); + + if (minimumOilPressure_.size() > 0) + minimumOilPressure_[globalDofIdx] = + std::min(Opm::getValue(fs.pressure(oilPhaseIdx)), + problem.minOilPressure(globalDofIdx)); + + if (overburdenPressure_.size() > 0) + overburdenPressure_[globalDofIdx] = problem.overburdenPressure(globalDofIdx); + + if (rockCompPorvMultiplier_.size() > 0) + rockCompPorvMultiplier_[globalDofIdx] = problem.template rockCompPoroMultiplier(intQuants, globalDofIdx); + + if (rockCompTransMultiplier_.size() > 0) + rockCompTransMultiplier_[globalDofIdx] = problem.template rockCompTransMultiplier(intQuants, globalDofIdx); + + const auto& matLawManager = problem.materialLawManager(); if (matLawManager->enableHysteresis()) { if (pcSwMdcOw_.size() > 0 && krnSwMdcOw_.size() > 0) { matLawManager->oilWaterHysteresisParams( @@ -525,7 +557,7 @@ public: // This can be removed when ebos has 100% controll over output if (elemCtx.simulator().episodeIndex() < 0 && FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) { - const auto& fsInitial = elemCtx.simulator().problem().initialFluidState(globalDofIdx); + const auto& fsInitial = problem.initialFluidState(globalDofIdx); // use initial rs and rv values if (rv_.size() > 0) @@ -833,6 +865,23 @@ public: if (bubblePointPressure_.size() > 0) sol.insert ("PBUB", Opm::UnitSystem::measure::pressure, std::move(bubblePointPressure_), Opm::data::TargetType::RESTART_AUXILIARY); + + if (swMax_.size() > 0) + sol.insert ("SWMAX", Opm::UnitSystem::measure::identity, std::move(swMax_), Opm::data::TargetType::RESTART_SOLUTION); + + if (minimumOilPressure_.size() > 0) + sol.insert ("PRESROCC", Opm::UnitSystem::measure::pressure, std::move(minimumOilPressure_), Opm::data::TargetType::RESTART_SOLUTION); + + if (overburdenPressure_.size() > 0) + sol.insert ("PRES_OVB", Opm::UnitSystem::measure::pressure, std::move(overburdenPressure_), Opm::data::TargetType::RESTART_SOLUTION); + + if (rockCompPorvMultiplier_.size() > 0) + sol.insert ("PORV_RC", Opm::UnitSystem::measure::identity, std::move(rockCompPorvMultiplier_), Opm::data::TargetType::RESTART_SOLUTION); + + if (rockCompTransMultiplier_.size() > 0) + sol.insert ("TMULT_RC", Opm::UnitSystem::measure::identity, std::move(rockCompTransMultiplier_), Opm::data::TargetType::RESTART_SOLUTION); + + // Fluid in place for (int i = 0; i 0) { @@ -1400,6 +1449,12 @@ private: ScalarBuffer ppcw_; ScalarBuffer bubblePointPressure_; ScalarBuffer dewPointPressure_; + ScalarBuffer rockCompPorvMultiplier_; + ScalarBuffer rockCompTransMultiplier_; + ScalarBuffer swMax_; + ScalarBuffer overburdenPressure_; + ScalarBuffer minimumOilPressure_; + std::vector failedCellsPb_; std::vector failedCellsPd_; std::vector fipnum_; diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 00b03445a..4aea50fdc 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -79,12 +79,17 @@ #include #include #include +#include +#include + #include #include #include #include #include #include +#include +#include #include #include @@ -420,6 +425,7 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper; typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; typedef typename GET_PROP_TYPE(TypeTag, EclWellModel) EclWellModel; typedef typename GET_PROP_TYPE(TypeTag, EclAquiferModel) EclAquiferModel; @@ -437,6 +443,8 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) typedef typename GridView::template Codim<0>::Iterator ElementIterator; + typedef Opm::UniformXTabulated2DFunction TabulatedTwoDFunction; + struct RockParams { Scalar referencePressure; Scalar compressibility; @@ -668,7 +676,7 @@ public: updatePffDofData_(); if (GET_PROP_VALUE(TypeTag, EnablePolymer)) { - const auto& vanguard = simulator.vanguard(); + const auto& vanguard = this->simulator().vanguard(); const auto& gridView = vanguard.gridView(); int numElements = gridView.size(/*codim=*/0); maxPolymerAdsorption_.resize(numElements, 0.0); @@ -763,7 +771,8 @@ public: // re-compute all quantities which may possibly be affected. transmissibilities_.update(); - updatePorosity_(); + referencePorosity_[1] = referencePorosity_[0]; + updateReferencePorosity_(); updatePffDofData_(); } @@ -825,6 +834,7 @@ public: { const auto& simulator = this->simulator(); int epsiodeIdx = simulator.episodeIndex(); + bool invalidateIntensiveQuantities = false; const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx); if (drsdtActive_()) // DRSDT is enabled @@ -834,7 +844,18 @@ public: if (drvdtActive_()) // DRVDT is enabled for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) - maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*simulator.timeStepSize(); + maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*this->simulator().timeStepSize(); + + if (enableExperiments) { + // update maximum water saturation and minimum pressure + // used when ROCKCOMP is activated + const bool invalidateFromMaxWaterSat = updateMaxWaterSaturation_(); + const bool invalidateFromMinPressure = updateMinPressure_(); + invalidateIntensiveQuantities = invalidateFromMaxWaterSat || invalidateFromMinPressure; + } + + if (invalidateIntensiveQuantities) + this->model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0); wellModel_.beginTimeStep(); aquiferModel_.beginTimeStep(); @@ -847,10 +868,11 @@ public: * term for the solution of the previous time step. * * For quite technical reasons, the storage term cannot be recycled if either DRSDT - * or DRVDT are active in ebos. + * or DRVDT are active in ebos. Nor if the porosity is changes between timesteps + * using a pore volume multiplier (i.e., poreVolumeMultiplier() != 1.0) */ bool recycleFirstIterationStorage() const - { return !drsdtActive_() && !drvdtActive_(); } + { return !drsdtActive_() && !drvdtActive_() && rockCompPoroMult_.empty(); } /*! * \brief Called by the simulator before each Newton-Raphson iteration. @@ -1061,23 +1083,30 @@ public: /*! * \copydoc FvBaseMultiPhaseProblem::porosity + * + * For the EclProblem, this method is identical to referencePorosity(). The intensive + * quantities object may apply various multipliers (e.g. ones which model rock + * compressibility and water induced rock compaction) to it which depend on the + * current physical conditions. */ template Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const { unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx); - return porosity_[globalSpaceIdx]; + return referencePorosity_[timeIdx][globalSpaceIdx]; } /*! * \brief Returns the porosity of an element * - * Note that this method is *not* part of the generic eWoms problem API because it - * would bake the assumption that only the elements are the degrees of freedom into - * the interface. + * The reference porosity of an element is the porosity of the medium before modified + * by the current solution. Note that this method is *not* part of the generic eWoms + * problem API because it would bake the assumption that only the elements are the + * degrees of freedom into the interface. */ - Scalar porosity(unsigned elementIdx) const - { return porosity_[elementIdx]; } + Scalar referencePorosity(unsigned elementIdx, unsigned timeIdx) const + { return referencePorosity_[timeIdx][elementIdx]; } + /*! * \brief Returns the depth of an degree of freedom [m] @@ -1550,8 +1579,10 @@ public: } /*! - * \brief Returns an element's maximum oil phase saturation observed during the - * simulation. + * \brief Returns an element's historic maximum oil phase saturation that was + * observed during the simulation. + * + * In this context, "historic" means the the time before the current timestep began. * * This is a bit of a hack from the conceptional point of view, but it is required to * match the results of the 'flow' and ECLIPSE 100 simulators. @@ -1568,6 +1599,8 @@ public: * \brief Sets an element's maximum oil phase saturation observed during the * simulation. * + * In this context, "historic" means the the time before the current timestep began. + * * This a hack on top of the maxOilSaturation() hack but it is currently required to * do restart externally. i.e. from the flow code. */ @@ -1579,6 +1612,43 @@ public: maxOilSaturation_[globalDofIdx] = value; } + + /*! + * \brief Returns an element's historic maximum water phase saturation that was + * observed during the simulation. + * + * In this context, "historic" means the the time before the current timestep began. + * + * This is used for output of the maximum water saturation used as input + * for water induced rock compation ROCK2D/ROCK2DTR. + */ + Scalar maxWaterSaturation(unsigned globalDofIdx) const + { + if (maxWaterSaturation_.empty()) + return 0.0; + + return maxWaterSaturation_[globalDofIdx]; + } + + + /*! + * \brief Returns an element's historic minimum pressure of the oil phase that was + * observed during the simulation. + * + * In this context, "historic" means the the time before the current timestep began. + * + * This is used for output of the minimum pressure used as input + * for the irreversible rock compation option. + */ + Scalar minOilPressure(unsigned globalDofIdx) const + { + if (minOilPressure_.empty()) + return 0.0; + + return minOilPressure_[globalDofIdx]; + } + + /*! * \brief Returns a reference to the ECL well manager used by the problem. * @@ -1673,6 +1743,83 @@ public: +std::to_string(double(simulator.timeStepSize()))); } + /*! + * \brief Calculate the porosity multiplier due to water induced rock compaction. + * + * TODO: The API of this is a bit ad-hoc, it would be better to use context objects. + */ + template + LhsEval rockCompPoroMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const + { + if (!enableExperiments || rockCompPoroMult_.size() == 0) + return 1.0; + + unsigned tableIdx = 0; + if (!rockTableIdx_.empty()) + 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 = + Opm::min(Opm::decay(fs.pressure(oilPhaseIdx)), + minOilPressure_[elementIdx]); + + if (!overburdenPressure_.empty()) + effectiveOilPressure -= overburdenPressure_[elementIdx]; + + return rockCompPoroMult_[tableIdx].eval(effectiveOilPressure, SwDeltaMax, /*extrapolation=*/true); + } + + /*! + * \brief Calculate the transmissibility multiplier due to water induced rock compaction. + * + * TODO: The API of this is a bit ad-hoc, it would be better to use context objects. + */ + template + LhsEval rockCompTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const + { + if (!enableExperiments || rockCompTransMult_.size() == 0) + return 1.0; + + unsigned tableIdx = 0; + if (!rockTableIdx_.empty()) + 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 = + Opm::min(Opm::decay(fs.pressure(oilPhaseIdx)), + minOilPressure_[elementIdx]); + + if (overburdenPressure_.size() > 0) + effectiveOilPressure -= overburdenPressure_[elementIdx]; + + return rockCompTransMult_[tableIdx].eval(effectiveOilPressure, SwDeltaMax, /*extrapolation=*/true); + } + + /*! + * \brief Get the pressure of the overburden. + * + * This method is mainly for output. + */ + Scalar overburdenPressure(unsigned elementIdx) const + { + if (!enableExperiments || overburdenPressure_.size() == 0) + return 0.0; + + return overburdenPressure_[elementIdx]; + } + private: void checkDeckCompatibility_() const @@ -1878,6 +2025,62 @@ private: return false; } + bool updateMaxWaterSaturation_() + { + // water compaction is activated in ROCKCOMP + if (maxWaterSaturation_.size()== 0) + return false; + + maxWaterSaturation_[/*timeIdx=*/1] = maxWaterSaturation_[/*timeIdx=*/0]; + ElementContext elemCtx(this->simulator()); + const auto& vanguard = this->simulator().vanguard(); + auto elemIt = vanguard.gridView().template begin(); + const auto& elemEndIt = vanguard.gridView().template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& elem = *elemIt; + + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + + unsigned compressedDofIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& iq = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& fs = iq.fluidState(); + + Scalar Sw = Opm::decay(fs.saturation(waterPhaseIdx)); + maxWaterSaturation_[compressedDofIdx] = std::max(maxWaterSaturation_[compressedDofIdx], Sw); + } + + return true; + } + + bool updateMinPressure_() + { + // IRREVERS option is used in ROCKCOMP + if (minOilPressure_.size() == 0) + return false; + + ElementContext elemCtx(this->simulator()); + const auto& vanguard = this->simulator().vanguard(); + auto elemIt = vanguard.gridView().template begin(); + const auto& elemEndIt = vanguard.gridView().template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& elem = *elemIt; + + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + + unsigned compressedDofIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& iq = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& fs = iq.fluidState(); + + minOilPressure_[compressedDofIdx] = + std::min(minOilPressure_[compressedDofIdx], + Opm::getValue(fs.pressure(oilPhaseIdx))); + } + + return true; + } + void readRockParameters_() { const auto& simulator = this->simulator(); @@ -1885,21 +2088,23 @@ private: const auto& eclState = simulator.vanguard().eclState(); const auto& vanguard = simulator.vanguard(); - // the ROCK keyword has not been specified, so we don't need - // to read rock parameters - if (!deck.hasKeyword("ROCK")) - return; - - const auto& rockKeyword = deck.getKeyword("ROCK"); - rockParams_.resize(rockKeyword.size()); - for (size_t rockRecordIdx = 0; rockRecordIdx < rockKeyword.size(); ++ rockRecordIdx) { - const auto& rockRecord = rockKeyword.getRecord(rockRecordIdx); - rockParams_[rockRecordIdx].referencePressure = - rockRecord.getItem("PREF").getSIDouble(0); - rockParams_[rockRecordIdx].compressibility = - rockRecord.getItem("COMPRESSIBILITY").getSIDouble(0); + // read the rock compressibility parameters + if (deck.hasKeyword("ROCK")) { + const auto& rockKeyword = deck.getKeyword("ROCK"); + rockParams_.resize(rockKeyword.size()); + for (size_t rockRecordIdx = 0; rockRecordIdx < rockKeyword.size(); ++ rockRecordIdx) { + const auto& rockRecord = rockKeyword.getRecord(rockRecordIdx); + rockParams_[rockRecordIdx].referencePressure = + rockRecord.getItem("PREF").getSIDouble(0); + rockParams_[rockRecordIdx].compressibility = + rockRecord.getItem("COMPRESSIBILITY").getSIDouble(0); + } } + // read the parameters for water-induced rock compaction + if (enableExperiments) + readRockCompactionParameters_(); + // check the kind of region which is supposed to be used by checking the ROCKOPTS // keyword. note that for some funny reason, the ROCK keyword uses PVTNUM by // default, *not* ROCKNUM! @@ -1920,20 +2125,144 @@ private: } } + // If ROCKCOMP is used and ROCKNUM is specified ROCK2D ROCK2DTR ROCKTAB etc. uses ROCKNUM + // to give the correct table index. + if (deck.hasKeyword("ROCKCOMP") && eclState.get3DProperties().hasDeckIntGridProperty("ROCKNUM")) + propName = "ROCKNUM"; + // the deck does not specify the selected keyword, so everything uses the first // record of ROCK. - if (!eclState.get3DProperties().hasDeckIntGridProperty(propName)) + if (eclState.get3DProperties().hasDeckIntGridProperty(propName)) { + const std::vector& tablenumData = + eclState.get3DProperties().getIntGridProperty(propName).getData(); + unsigned numElem = vanguard.gridView().size(0); + rockTableIdx_.resize(numElem); + for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) { + unsigned cartElemIdx = vanguard.cartesianIndex(elemIdx); + + // reminder: Eclipse uses FORTRAN-style indices + rockTableIdx_[elemIdx] = tablenumData[cartElemIdx] - 1; + } + } + + // Store overburden pressure pr element + const auto& overburdTables = eclState.getTableManager().getOverburdTables(); + if (!overburdTables.empty()) { + unsigned numElem = vanguard.gridView().size(0); + overburdenPressure_.resize(numElem,0.0); + + const auto& rockcomp = deck.getKeyword("ROCKCOMP"); + const auto& rockcompRecord = rockcomp.getRecord(0); + size_t numRocktabTables = rockcompRecord.getItem("NTROCC").template get< int >(0); + + if (overburdTables.size() != numRocktabTables) + throw std::runtime_error(std::to_string(numRocktabTables) +" OVERBURD tables is expected, but " + std::to_string(overburdTables.size()) +" is provided"); + + std::vector> overburdenTables(numRocktabTables); + for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { + const Opm::OverburdTable& overburdTable = overburdTables.template getTable(regionIdx); + overburdenTables[regionIdx].setXYContainers(overburdTable.getDepthColumn(),overburdTable.getOverburdenPressureColumn()); + } + + for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) { + unsigned tableIdx = 0; + if (!rockTableIdx_.empty()) { + tableIdx = rockTableIdx_[elemIdx]; + } + overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(elementCenterDepth_[elemIdx], /*extrapolation=*/true); + } + } + + } + + void readRockCompactionParameters_() + { + const auto& vanguard = this->simulator().vanguard(); + const auto& deck = vanguard.deck(); + const auto& eclState = vanguard.eclState(); + + if (!deck.hasKeyword("ROCKCOMP")) + return; // deck does not enable rock compaction + + const auto& rockcomp = deck.getKeyword("ROCKCOMP"); + //for (size_t rockRecordIdx = 0; rockRecordIdx < rockcomp.size(); ++ rockRecordIdx) { + assert(rockcomp.size() == 1); + const auto& rockcompRecord = rockcomp.getRecord(0); + const auto& option = rockcompRecord.getItem("HYSTERESIS").getTrimmedString(0); + if (option == "REVERS") { + // interpolate the porv volume multiplier using the pressure in the cell + } + else if (option == "IRREVERS") { + // interpolate the porv volume multiplier using the minimum pressure in the cell + // i.e. don't allow re-inflation. + unsigned numElem = vanguard.gridView().size(0); + minOilPressure_.resize(numElem, 1e99); + } + else if (option == "NO") + // rock compaction turned on but disabled by ROCKCOMP option return; + else + throw std::runtime_error("ROCKCOMP option " + option + " not supported for item 1"); - const std::vector& tablenumData = - eclState.get3DProperties().getIntGridProperty(propName).getData(); - unsigned numElem = vanguard.gridView().size(0); - rockTableIdx_.resize(numElem); - for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) { - unsigned cartElemIdx = vanguard.cartesianIndex(elemIdx); + size_t numRocktabTables = rockcompRecord.getItem("NTROCC").template get(0); + const auto& waterCompactionItem = rockcompRecord.getItem("WATER_COMPACTION").getTrimmedString(0); + bool waterCompaction = false; + if (waterCompactionItem == "YES") { + waterCompaction = true; + unsigned numElem = vanguard.gridView().size(0); + maxWaterSaturation_.resize(numElem, 0.0); + } + else + throw std::runtime_error("ROCKCOMP option " + waterCompactionItem + " not supported for item 3. Only YES is supported"); - // reminder: Eclipse uses FORTRAN-style indices - rockTableIdx_[elemIdx] = tablenumData[cartElemIdx] - 1; + if (waterCompaction) { + const auto& rock2dTables = eclState.getTableManager().getRock2dTables(); + const auto& rock2dtrTables = eclState.getTableManager().getRock2dtrTables(); + const auto& rockwnodTables = eclState.getTableManager().getRockwnodTables(); + + if (rock2dTables.size() != numRocktabTables) + throw std::runtime_error("Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables) + +" ROCK2D tables is expected, but " + std::to_string(rock2dTables.size()) +" is provided"); + + if (rockwnodTables.size() != numRocktabTables) + 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)); + for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) { + const Opm::RockwnodTable& rockwnodTable = rockwnodTables.template getTable(regionIdx); + const auto& rock2dTable = rock2dTables[regionIdx]; + + if (rockwnodTable.getSaturationColumn().size() != rock2dTable.sizeMultValues()) + 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)); + for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx) + rockCompPoroMult_[regionIdx].appendSamplePoint(xIdx, + rockwnodTable.getSaturationColumn()[yIdx], + rock2dTable.getPvmultValue(xIdx, yIdx)); + } + } + + if (!rock2dtrTables.empty()) { + rockCompTransMult_.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]; + + if (rockwnodTable.getSaturationColumn().size() != rock2dtrTable.sizeMultValues()) + 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)); + for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx) + rockCompTransMult_[regionIdx].appendSamplePoint(xIdx, + rockwnodTable.getSaturationColumn()[yIdx], + rock2dtrTable.getTransMultValue(xIdx, yIdx)); + } + } + } } } @@ -1955,7 +2284,8 @@ private: //////////////////////////////// // porosity - updatePorosity_(); + updateReferencePorosity_(); + referencePorosity_[1] = referencePorosity_[0]; //////////////////////////////// //////////////////////////////// @@ -1990,7 +2320,7 @@ private: thermalLawManager_->initFromDeck(deck, eclState, compressedToCartesianElemIdx); } - void updatePorosity_() + void updateReferencePorosity_() { const auto& simulator = this->simulator(); const auto& vanguard = simulator.vanguard(); @@ -2000,7 +2330,7 @@ private: size_t numDof = this->model().numGridDof(); - porosity_.resize(numDof); + referencePorosity_[/*timeIdx=*/0].resize(numDof); const std::vector& porvData = props.getDoubleGridProperty("PORV").getData(); @@ -2043,7 +2373,7 @@ private: // be larger than 1.0! Scalar dofVolume = simulator.model().dofTotalVolume(dofIdx); assert(dofVolume > 0.0); - porosity_[dofIdx] = poreVolume/dofVolume; + referencePorosity_[/*timeIdx=*/0][dofIdx] = poreVolume/dofVolume; } } @@ -2069,6 +2399,19 @@ private: readBlackoilExtentionsInitialConditions_(); + //initialize min/max values + size_t numElems = this->model().numGridDof(); + for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) { + const auto& fs = initialFluidStates_[elemIdx]; + if(maxWaterSaturation_.size() > 0) + maxWaterSaturation_[elemIdx] = std::max(maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx)); + if(maxOilSaturation_.size() > 0) + maxOilSaturation_[elemIdx] = std::max(maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx)); + if(minOilPressure_.size() > 0) + minOilPressure_[elemIdx] = std::min(minOilPressure_[elemIdx], fs.pressure(oilPhaseIdx)); + } + + } void readEquilInitialCondition_() @@ -2712,6 +3055,10 @@ private: // first thing in the morning, limit the time step size to the maximum size dtNext = std::min(dtNext, maxTimeStepSize_); + // use at least slightly more than half of the maximum time step size by default + if (dtNext < maxTimeStepSize_ && maxTimeStepSize_ < dtNext*2) + dtNext = 1.01 * maxTimeStepSize_/2.0; + Scalar remainingEpisodeTime = simulator.episodeStartTime() + simulator.episodeLength() - (simulator.startTime() + simulator.time()); @@ -2740,7 +3087,7 @@ private: static std::string briefDescription_; - std::vector porosity_; + std::array, 2> referencePorosity_; std::vector elementCenterDepth_; EclTransmissibility transmissibilities_; @@ -2773,6 +3120,12 @@ private: std::vector maxDRv_; constexpr static Scalar freeGasMinSaturation_ = 1e-7; std::vector maxOilSaturation_; + std::vector maxWaterSaturation_; + std::vector overburdenPressure_; + std::vector minOilPressure_; + + std::vector rockCompPoroMult_; + std::vector rockCompTransMult_; bool enableDriftCompensation_; GlobalEqVector drift_; diff --git a/ebos/eclthresholdpressure.hh b/ebos/eclthresholdpressure.hh index 7b17aee86..d972e70a6 100644 --- a/ebos/eclthresholdpressure.hh +++ b/ebos/eclthresholdpressure.hh @@ -252,9 +252,9 @@ private: continue; // don't include connections with negligible flow - const Scalar& trans = simulator_.problem().transmissibility(elemCtx, i, j); - const Scalar& faceArea = face.area(); - if ( std::abs(faceArea * trans) < 1e-18) + const Evaluation& trans = simulator_.problem().transmissibility(elemCtx, i, j); + Scalar faceArea = face.area(); + if (std::abs(faceArea*Opm::getValue(trans)) < 1e-18) continue; // determine the maximum difference of the pressure of any phase over the diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 0ec671752..28ecec9d2 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -591,7 +591,7 @@ namespace Opm { const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& fs = intQuants.fluidState(); - const double pvValue = ebosProblem.porosity(cell_idx) * ebosModel.dofTotalVolume( cell_idx ); + const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx ); pvSumLocal += pvValue; for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 44cb377ca..1625e59d8 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -319,6 +319,7 @@ namespace Opm void computePerfRate(const IntensiveQuantities& intQuants, const std::vector& mob, const EvalWell& bhp, + const double Tw, const int perf, const bool allow_cf, std::vector& cq_s, diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 9b5e60d1d..3bf520a0c 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -282,6 +282,7 @@ namespace Opm computePerfRate(const IntensiveQuantities& intQuants, const std::vector& mob, const EvalWell& bhp, + const double Tw, const int perf, const bool allow_cf, std::vector& cq_s, @@ -317,7 +318,6 @@ namespace Opm return; } - const double Tw = well_index_[perf]; // compute component volumetric rates at standard conditions for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { const EvalWell cq_p = - Tw * (mob[componentIdx] * drawdown); @@ -355,7 +355,6 @@ namespace Opm } // injection perforations total volume rates - const double Tw = well_index_[perf]; const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown); // surface volume fraction of fluids within wellbore @@ -487,7 +486,9 @@ namespace Opm std::vector cq_s(num_components_, 0.0); double perf_dis_gas_rate = 0.; double perf_vap_oil_rate = 0.; - computePerfRate(intQuants, mob, bhp, perf, allow_cf, + double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); + const double Tw = well_index_[perf] * trans_mult; + computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf, cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); // updating the solution gas rate and solution oil rate @@ -1274,7 +1275,7 @@ namespace Opm } // the well index associated with the connection - const double tw_perf = well_index_[perf]; + const double tw_perf = well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier(int_quantities, cell_idx); // TODO: there might be some indices related problems here // phases vs components @@ -2203,11 +2204,13 @@ namespace Opm // flux for each perforation std::vector mob(num_components_, 0.0); getMobility(ebosSimulator, perf, mob, deferred_logger); + double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); + const double Tw = well_index_[perf] * trans_mult; std::vector cq_s(num_components_, 0.0); double perf_dis_gas_rate = 0.; double perf_vap_oil_rate = 0.; - computePerfRate(intQuants, mob, bhp, perf, allow_cf, + computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf, cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); for(int p = 0; p < np; ++p) { @@ -2610,7 +2613,9 @@ namespace Opm std::vector cq_s(num_components_,0.0); double perf_dis_gas_rate = 0.; double perf_vap_oil_rate = 0.; - computePerfRate(int_quant, mob, bhp, perf, allow_cf, + double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier(int_quant, cell_idx); + const double Tw = well_index_[perf] * trans_mult; + computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); // TODO: make area a member const double area = 2 * M_PI * perf_rep_radius_[perf] * perf_length_[perf];