Add support for drsdtr and drvdtr

This PR also adds possibility for schedule dependent drsdt values

ebos/eclbasevanguard.hh

@@ -240,7 +240,7 @@ public:
                 const auto& grid = internalEclState_->getInputGrid();
                 const Opm::TableManager table ( *internalDeck_ );
                 const Opm::Eclipse3DProperties eclipseProperties (*internalDeck_  , table, grid);
-                internalSchedule_.reset(new Opm::Schedule(*internalDeck_, grid, eclipseProperties, Opm::Phases(true, true, true), parseContext ));
+                internalSchedule_.reset(new Opm::Schedule(*internalDeck_, grid, eclipseProperties, internalEclState_->runspec(), parseContext ));
                 internalSummaryConfig_.reset(new Opm::SummaryConfig(*internalDeck_, *internalSchedule_, table,  parseContext));
             }
 

ebos/eclproblem.hh

@@ -523,72 +523,17 @@ public:
         // or DRVDT up to the current time step in the schedule section counts, presence
         // of VAPPARS alone is not sufficient to disable DR[SV]DT. TODO: implment support
         // for this in opm-parser's Schedule object"
-        drsdtActive_ = false;
-        drvdtActive_ = false;
-        vapparsActive_ = false;
-
-        if (deck.hasKeyword("VAPPARS")) {
-            vapparsActive_ = true;
-
-            size_t numDof = this->model().numGridDof();
-            maxOilSaturation_.resize(numDof, 0.0);
-
-            // TODO: update the PVT objects. this is only required if VAPPARS becomes a
-            // fully dynamic keyword.
-        }
 
         // deal with DRSDT
         const auto& eclState = simulator.vanguard().eclState();
         unsigned ntpvt = eclState.runspec().tabdims().getNumPVTTables();
-        maxDRsDt_.resize(ntpvt, 0.0);
+        maxDRs_.resize(ntpvt, 1e30);
-        maxDRs_.resize(ntpvt, -1.0);
         dRsDtOnlyFreeGas_.resize(ntpvt, false);
-        if (!vapparsActive_) {
+        size_t numDof = this->model().numGridDof();
-            if (deck.hasKeyword("DRSDT")) {
+        lastRs_.resize(numDof, 0.0);
-                drsdtActive_ = !vapparsActive_;
+        maxDRv_.resize(ntpvt, 1e30);
-                const auto& drsdtKeyword = deck.getKeyword("DRSDT");
+        lastRv_.resize(numDof, 0.0);
-                std::fill(maxDRsDt_.begin(), maxDRsDt_.end(), drsdtKeyword.getRecord(0).getItem("DRSDT_MAX").getSIDouble(0));
+        maxOilSaturation_.resize(numDof, 0.0);
-                size_t numDof = this->model().numGridDof();
-                lastRs_.resize(numDof, 0.0);
-                std::string drsdtFlag =
-                        drsdtKeyword.getRecord(0).getItem("Option").getTrimmedString(0);
-                std::transform(drsdtFlag.begin(), drsdtFlag.end(), drsdtFlag.begin(), ::toupper);
-
-                std::fill(dRsDtOnlyFreeGas_.begin(), dRsDtOnlyFreeGas_.end(), (drsdtFlag == "FREE"));
-            } else if (deck.hasKeyword("DRSDTR")) {
-                drsdtActive_ = !vapparsActive_;
-                const auto& drsdtrKeyword = deck.getKeyword("DRSDTR");
-                size_t numDof = this->model().numGridDof();
-                lastRs_.resize(numDof, 0.0);
-                for (size_t recordIdx = 0; recordIdx < drsdtrKeyword.size(); ++recordIdx) {
-                    const Opm::DeckRecord& record = drsdtrKeyword.getRecord(recordIdx);
-                    maxDRsDt_[recordIdx] = record.getItem("DRSDT_MAX").getSIDouble(0);
-                    std::string drsdtFlag = record.getItem("Option").getTrimmedString(0);
-                    std::transform(drsdtFlag.begin(), drsdtFlag.end(), drsdtFlag.begin(), ::toupper);
-                    dRsDtOnlyFreeGas_[recordIdx] = (drsdtFlag == "FREE");
-                }
-            }
-        }
-
-        // deal with DRVDT
-        maxDRvDt_.resize(ntpvt, 0.0);
-        maxDRv_.resize(ntpvt, -1.0);
-        if (!vapparsActive_) {
-            if (deck.hasKeyword("DRVDT")) {
-                const auto& drvdtKeyword = deck.getKeyword("DRVDT");
-                std::fill(maxDRvDt_.begin(), maxDRvDt_.end(), drvdtKeyword.getRecord(0).getItem("DRVDT_MAX").getSIDouble(0));
-                size_t numDof = this->model().numGridDof();
-                lastRv_.resize(numDof, 0.0);
-            } else if (deck.hasKeyword("DRVDTR")) {
-                const auto& drvdtrKeyword = deck.getKeyword("DRVDTR");
-                size_t numDof = this->model().numGridDof();
-                lastRv_.resize(numDof, 0.0);
-                for (size_t recordIdx = 0; recordIdx < drvdtrKeyword.size(); ++recordIdx) {
-                    const Opm::DeckRecord& record = drvdtrKeyword.getRecord(recordIdx);
-                    maxDRvDt_[recordIdx] = record.getItem("DRVDT_MAX").getSIDouble(0);
-                }
-            }
-        }
 
         initFluidSystem_();
         updateElementDepths_();
@@ -760,18 +705,21 @@ public:
      */
     void beginTimeStep()
     {
-        if (drsdtActive_)
+        int epsiodeIdx = this->simulator().episodeIndex();
+        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+        if (drsdtActive_())
             // DRSDT is enabled
             for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx )
-                maxDRs_[pvtRegionIdx] = maxDRsDt_[pvtRegionIdx]*this->simulator().timeStepSize();
+                maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*this->simulator().timeStepSize();
 
-        if (drvdtActive_)
+        if (drvdtActive_())
             // DRVDT is enabled
             for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx )
-                maxDRv_[pvtRegionIdx] = maxDRvDt_[pvtRegionIdx]*this->simulator().timeStepSize();
+                maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*this->simulator().timeStepSize();
 
-        if (!GET_PROP_VALUE(TypeTag, DisableWells))
+        if (!GET_PROP_VALUE(TypeTag, DisableWells)) {
             wellModel_.beginTimeStep();
+        }
 
         aquiferModel_.beginTimeStep();
     }
@@ -834,6 +782,7 @@ public:
             initialFluidStates_.clear();
         }
 
+
         updateCompositionChangeLimits_();
     }
 
@@ -853,11 +802,6 @@ public:
             simulator.setFinished(true);
             return;
         }
-
-        if (!GET_PROP_VALUE(TypeTag, DisableWells))
-            wellModel_.endEpisode();
-
-        aquiferModel_.endEpisode();
     }
 
     /*!
@@ -1325,6 +1269,7 @@ public:
      */
     void initialSolutionApplied()
     {
+
         if (!GET_PROP_VALUE(TypeTag, DisableWells)) {
             // initialize the wells. Note that this needs to be done after initializing the
             // intrinsic permeabilities and the after applying the initial solution because
@@ -1382,7 +1327,7 @@ public:
     Scalar maxGasDissolutionFactor(unsigned globalDofIdx) const
     {
         int pvtRegionIdx = pvtRegionIndex(globalDofIdx);
-        if (!drsdtActive_ || maxDRs_[pvtRegionIdx] < 0.0)
+        if (!drsdtActive_() || maxDRs_[pvtRegionIdx] < 0.0)
             return std::numeric_limits<Scalar>::max()/2;
 
         return lastRs_[globalDofIdx] + maxDRs_[pvtRegionIdx];
@@ -1395,7 +1340,7 @@ public:
     Scalar maxOilVaporizationFactor(unsigned globalDofIdx) const
     {
         int pvtRegionIdx = pvtRegionIndex(globalDofIdx);
-        if (!drvdtActive_ || maxDRv_[pvtRegionIdx] < 0.0)
+        if (!drvdtActive_() || maxDRv_[pvtRegionIdx] < 0.0)
             return std::numeric_limits<Scalar>::max()/2;
 
         return lastRv_[globalDofIdx] + maxDRv_[pvtRegionIdx];
@@ -1410,7 +1355,7 @@ public:
      */
     Scalar maxOilSaturation(unsigned globalDofIdx) const
     {
-        if (!vapparsActive_)
+        if (!vapparsActive())
             return 0.0;
 
         return maxOilSaturation_[globalDofIdx];
@@ -1425,7 +1370,7 @@ public:
      */
     void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
     {
-        if (!vapparsActive_)
+        if (!vapparsActive())
             return;
 
         maxOilSaturation_[globalDofIdx] = value;
@@ -1451,10 +1396,26 @@ public:
 
     bool vapparsActive() const
     {
-        return vapparsActive_;
+        int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0 );
+        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+        return (oilVaporizationControl.getType() == Opm::OilVaporizationEnum::VAPPARS);
     }
 
 private:
+    bool drsdtActive_() const
+    {
+        int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0 );
+        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+        return (oilVaporizationControl.drsdtActive());
+
+    }
+    bool drvdtActive_() const
+    {
+        int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0 );
+        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+        return (oilVaporizationControl.drvdtActive());
+
+    }
     Scalar cellCenterDepth( const Element& element ) const
     {
         typedef typename Element :: Geometry Geometry;
@@ -1495,7 +1456,10 @@ private:
     {
         // update the "last Rs" values for all elements, including the ones in the ghost
         // and overlap regions
-        if (drsdtActive_) {
+        int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0 );
+        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+
+        if (oilVaporizationControl.drsdtActive()) {
             ElementContext elemCtx(this->simulator());
             const auto& vanguard = this->simulator().vanguard();
             auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
@@ -1513,7 +1477,7 @@ private:
                 typedef typename std::decay<decltype(fs) >::type FluidState;
 
                 int pvtRegionIdx = pvtRegionIndex(compressedDofIdx);
-                if (!dRsDtOnlyFreeGas_[pvtRegionIdx] || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_)
+                if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_)
                     lastRs_[compressedDofIdx] =
                         Opm::BlackOil::template getRs_<FluidSystem,
                                                        FluidState,
@@ -1525,7 +1489,7 @@ private:
 
         // update the "last Rv" values for all elements, including the ones in the ghost
         // and overlap regions
-        if (drvdtActive_) {
+        if (drvdtActive_()) {
             ElementContext elemCtx(this->simulator());
             const auto& vanguard = this->simulator().vanguard();
             auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
@@ -1553,7 +1517,7 @@ private:
     bool updateMaxOilSaturation_()
     {
         // we use VAPPARS
-        if (vapparsActive_) {
+        if (vapparsActive()) {
             ElementContext elemCtx(this->simulator());
             const auto& vanguard = this->simulator().vanguard();
             auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
@@ -1973,6 +1937,8 @@ private:
         }
     }
 
+
+
     // update the hysteresis parameters of the material laws for the whole grid
     bool updateHysteresis_()
     {
@@ -2152,18 +2118,12 @@ private:
     std::vector<Scalar> polymerConcentration_;
     std::vector<Scalar> solventSaturation_;
 
-    bool drsdtActive_; // if no, VAPPARS *might* be active
     std::vector<bool> dRsDtOnlyFreeGas_; // apply the DRSDT rate limit only to cells that exhibit free gas
     std::vector<Scalar> lastRs_;
-    std::vector<Scalar> maxDRsDt_;
     std::vector<Scalar> maxDRs_;
-    bool drvdtActive_; // if no, VAPPARS *might* be active
     std::vector<Scalar> lastRv_;
-    std::vector<Scalar> maxDRvDt_;
     std::vector<Scalar> maxDRv_;
     constexpr static Scalar freeGasMinSaturation_ = 1e-7;
-
-    bool vapparsActive_; // if no, DRSDT and/or DRVDT *might* be active
     std::vector<Scalar> maxOilSaturation_;
 
     EclWellModel wellModel_;