move the maximum oil saturation stuff into the problem

this also simplifies creating a more proper selection logic for the
DR[SV]DT keywords or VAPPARS.

ebos/ecloutputblackoilmodule.hh

@@ -271,14 +271,14 @@ public:
 
             }
             if (gasDissolutionFactor_.size() > 0) {
-                Scalar SoMax = elemCtx.model().maxOilSaturation(globalDofIdx);
+                Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx);
                 gasDissolutionFactor_[globalDofIdx] =
                     FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx, SoMax);
                 Opm::Valgrind::CheckDefined(gasDissolutionFactor_[globalDofIdx]);
 
             }
             if (oilVaporizationFactor_.size() > 0) {
-                Scalar SoMax = elemCtx.model().maxOilSaturation(globalDofIdx);
+                Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx);
                 oilVaporizationFactor_[globalDofIdx] =
                     FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(fs, gasPhaseIdx, pvtRegionIdx, SoMax);
                 Opm::Valgrind::CheckDefined(oilVaporizationFactor_[globalDofIdx]);
@@ -375,7 +375,7 @@ public:
             }
 
             if (soMax_.size() > 0)
-                soMax_[globalDofIdx] = elemCtx.simulator().model().maxOilSaturation(globalDofIdx);
+                soMax_[globalDofIdx] = elemCtx.simulator().problem().maxOilSaturation(globalDofIdx);
 
             const auto& matLawManager = elemCtx.simulator().problem().materialLawManager();
             if (matLawManager->enableHysteresis()) {
@@ -723,7 +723,7 @@ public:
     void initHysteresisParams(Simulator& simulator, unsigned elemIdx) const {
 
         if (soMax_.size() > 0)
-            simulator.model().setMaxOilSaturation(soMax_[elemIdx], elemIdx);
+            simulator.problem().setMaxOilSaturation(soMax_[elemIdx], elemIdx);
 
         if (simulator.problem().materialLawManager()->enableHysteresis()) {
             auto matLawManager = simulator.problem().materialLawManager();

ebos/eclproblem.hh

@@ -362,12 +362,25 @@ 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"
-        bool hasVappars = deck.hasKeyword("VAPPARS");
+        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
         maxDRsDt_ = 0.0;
         maxDRs_ = -1.0;
-        if (!hasVappars && deck.hasKeyword("DRSDT")) {
+        if (!vapparsActive_ && deck.hasKeyword("DRSDT")) {
+            drsdtActive_ = !vapparsActive_;
             const auto& drsdtKeyword = deck.getKeyword("DRSDT");
             maxDRsDt_ = drsdtKeyword.getRecord(0).getItem("DRSDT_MAX").getSIDouble(0);
             size_t numDof = this->model().numGridDof();
@@ -377,7 +390,7 @@ public:
         // deal with DRVDT
         maxDRvDt_ = 0.0;
         maxDRv_ = -1.0;
-        if (!hasVappars && deck.hasKeyword("DRVDT")) {
+        if (!vapparsActive_ && deck.hasKeyword("DRVDT")) {
             const auto& drvdtKeyword = deck.getKeyword("DVSDT");
             maxDRvDt_ = drvdtKeyword.getRecord(0).getItem("DRVDT_MAX").getSIDouble(0);
             size_t numDof = this->model().numGridDof();
@@ -502,10 +515,8 @@ public:
             simulator.setTimeStepSize(dt);
         }
 
-        if (updateHysteresis_())
+        bool doInvalidate = updateHysteresis_();
-            this->model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
+        doInvalidate = updateMaxOilSaturation_();
-
-        this->model().updateMaxOilSaturations();
 
         if (GET_PROP_VALUE(TypeTag, EnablePolymer))
             updateMaxPolymerAdsorption_();
@@ -514,6 +525,9 @@ public:
             // set up the wells
             wellManager_.beginEpisode(this->simulator().gridManager().eclState(),
                                       this->simulator().gridManager().schedule(), isOnRestart);
+
+        if (doInvalidate)
+            this->model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
     }
 
     /*!
@@ -1054,9 +1068,9 @@ public:
      * \brief Returns the maximum value of the gas dissolution factor at the current time
      *        for a given degree of freedom.
      */
-    Scalar maxGasDissolutionFactor(unsigned globalDofIdx OPM_UNUSED) const
+    Scalar maxGasDissolutionFactor(unsigned globalDofIdx) const
     {
-        if (lastRs_.empty() || maxDRs_ < 0.0)
+        if (!drsdtActive_ || maxDRs_ < 0.0)
             return std::numeric_limits<Scalar>::max()/2;
 
         return lastRs_[globalDofIdx] + maxDRs_;
@@ -1066,14 +1080,46 @@ public:
      * \brief Returns the maximum value of the oil vaporization factor at the current
      *        time for a given degree of freedom.
      */
-    Scalar maxOilVaporizationFactor(unsigned globalDofIdx OPM_UNUSED) const
+    Scalar maxOilVaporizationFactor(unsigned globalDofIdx) const
     {
-        if (lastRv_.empty() || maxDRv_ < 0.0)
+        if (!drvdtActive_ || maxDRv_ < 0.0)
             return std::numeric_limits<Scalar>::max()/2;
 
         return lastRv_[globalDofIdx] + maxDRv_;
     }
 
+    /*!
+     * \brief Returns an element's maximum oil phase saturation observed during the
+     *        simulation.
+     *
+     * 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.
+     */
+    Scalar maxOilSaturation(unsigned globalDofIdx) const
+    {
+        if (!vapparsActive_)
+            return 0.0;
+
+        return maxOilSaturation_[globalDofIdx];
+    }
+
+    /*!
+     * \brief Sets an element's maximum oil phase saturation observed during the
+     *        simulation.
+     *
+     * This a hack on top of the maxOilSaturation() hack but it is currently required to
+     * do restart externally. i.e. from the flow code.
+     *
+     * TODO: move the restart-from-ECL-restart-files functionality to EclProblem!
+     */
+    void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
+    {
+        if (!vapparsActive_)
+            return;
+
+        maxOilSaturation_[globalDofIdx] = value;
+    }
+
     /*!
      * \brief Returns a reference to the ECL well manager used by the problem.
      *
@@ -1135,7 +1181,7 @@ private:
     {
         // update the "last Rs" values for all elements, including the ones in the ghost
         // and overlap regions
-        if (!lastRs_.empty()) {
+        if (drsdtActive_) {
             ElementContext elemCtx(this->simulator());
             const auto& gridManager = this->simulator().gridManager();
             auto elemIt = gridManager.gridView().template begin</*codim=*/0>();
@@ -1162,7 +1208,7 @@ private:
 
         // update the "last Rv" values for all elements, including the ones in the ghost
         // and overlap regions
-        if (!lastRv_.empty()) {
+        if (drvdtActive_) {
             ElementContext elemCtx(this->simulator());
             const auto& gridManager = this->simulator().gridManager();
             auto elemIt = gridManager.gridView().template begin</*codim=*/0>();
@@ -1186,12 +1232,37 @@ private:
                                                                iq.pvtRegionIndex());
             }
         }
+    }
 
-        if (!lastRs_.empty() || !lastRv_.empty()) {
+    bool updateMaxOilSaturation_()
-            // we need to invalidate the intensive quantities cache here because the Rs
+    {
-            // values of the intensive quantities are potentially different now.
+        // we use VAPPARS
-            this->model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
+        if (vapparsActive_) {
+            ElementContext elemCtx(this->simulator());
+            const auto& gridManager = this->simulator().gridManager();
+            auto elemIt = gridManager.gridView().template begin</*codim=*/0>();
+            const auto& elemEndIt = gridManager.gridView().template end</*codim=*/0>();
+            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 So = Opm::decay<Scalar>(fs.saturation(oilPhaseIdx));
+
+                maxOilSaturation_[compressedDofIdx] = std::max(maxOilSaturation_[compressedDofIdx], So);
+            }
+
+            // we need to invalidate the intensive quantities cache here because the
+            // derivatives of Rs and Rv will most likely have changed
+            return true;
         }
+
+        return false;
     }
 
     void readRockParameters_()
@@ -1542,6 +1613,7 @@ private:
 
         }
     }
+
     void readBlackoilExtentionsInitialConditions_()
     {
         const auto& gridManager = this->simulator().gridManager();
@@ -1744,14 +1816,19 @@ private:
     std::vector<Scalar> polymerConcentration_;
     std::vector<Scalar> solventSaturation_;
 
+    bool drsdtActive_; // if no, VAPPARS *might* be active
     std::vector<Scalar> lastRs_;
     Scalar maxDRsDt_;
     Scalar maxDRs_;
 
+    bool drvdtActive_; // if no, VAPPARS *might* be active
     std::vector<Scalar> lastRv_;
     Scalar maxDRvDt_;
     Scalar maxDRv_;
 
+    bool vapparsActive_; // if no, DRSDT and/or DRVDT *might* be active
+    std::vector<Scalar> maxOilSaturation_;
+
     EclWellManager<TypeTag> wellManager_;
 
     std::unique_ptr< EclWriterType > eclWriter_;