Merge pull request #5527 from totto82/pvtnum_drdsdtcon

Support deactivation of individual DRSDT/DRVDT/DRSDTCON regions

opm/simulators/flow/FlowProblem.hpp

@@ -546,6 +546,10 @@ public:
         }
 
         bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex());
+        this->mixControls_.init(this->model().numGridDof(),
+                                episodeIdx,
+                                eclState.runspec().tabdims().getNumPVTTables());
+
 
         // set up the wells for the next episode.
         wellModel_.beginEpisode();
@@ -1943,8 +1947,7 @@ protected:
                                                             this->gravity_[dim - 1],
                                                             perm[dim - 1][dim - 1],
                                                             distZ,
-                                                            pvtRegionIdx,
+                                                            pvtRegionIdx);
-                                                            active);
                               }
             );
 

opm/simulators/flow/MixingRateControls.cpp

@@ -97,20 +97,19 @@ template<class FluidSystem>
 void MixingRateControls<FluidSystem>::
 init(std::size_t numDof, int episodeIdx, const unsigned ntpvt)
 {
-    // deal with DRSDT
+    // allocate DRSDT related vectors
-    //TODO We may want to only allocate these properties only if active.
+    if (this->drsdtActive(episodeIdx) && maxDRs_.empty()) {
-    //But since they may be activated at later time we need some more
-    //intrastructure to handle it
-
-        maxDRv_.resize(ntpvt, 1e30);
-        lastRv_.resize(numDof, 0.0);
         maxDRs_.resize(ntpvt, 1e30);
         dRsDtOnlyFreeGas_.resize(ntpvt, false);
         lastRs_.resize(numDof, 0.0);
+    }
+    if (this->drvdtActive(episodeIdx && maxDRv_.empty())) {
+        lastRv_.resize(numDof, 0.0);
         maxDRv_.resize(ntpvt, 1e30);
-        if (this->drsdtConvective(episodeIdx)) {
+    }
-            convectiveDrs_.resize(numDof, 1.0);
+    if (this->drsdtConvective(episodeIdx) && convectiveDrs_.empty()) {
-        }
+        convectiveDrs_.resize(numDof, 1.0);
+    }
 }
 
 template<class FluidSystem>
@@ -137,22 +136,49 @@ drsdtConvective(int episodeIdx) const
     return (oilVaporizationControl.drsdtConvective());
 }
 
+template<class FluidSystem>
+bool MixingRateControls<FluidSystem>::
+drsdtActive(int episodeIdx, std::size_t pvtRegionIdx) const
+{
+    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+    return (oilVaporizationControl.drsdtActive(pvtRegionIdx));
+}
+
+template<class FluidSystem>
+bool MixingRateControls<FluidSystem>::
+drvdtActive(int episodeIdx, std::size_t pvtRegionIdx) const
+{
+    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+    return (oilVaporizationControl.drvdtActive(pvtRegionIdx));
+}
+
+template<class FluidSystem>
+bool MixingRateControls<FluidSystem>::
+drsdtConvective(int episodeIdx, std::size_t pvtRegionIdx) const
+{
+    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+    return (oilVaporizationControl.drsdtConvective(pvtRegionIdx));
+}
+
+
 template<class FluidSystem>
 void MixingRateControls<FluidSystem>::
 updateExplicitQuantities(const int episodeIdx,
                          const Scalar timeStepSize)
 {
     const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    if (this->drsdtActive(episodeIdx)) {
+    // DRSDT is enabled
-        // DRSDT is enabled
+    for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) {
-        for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx)
+        if (oilVaporizationControl.drsdtActive(pvtRegionIdx)) {
             maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx) * timeStepSize;
+        }
     }
 
-    if (this->drvdtActive(episodeIdx)) {
+    // DRVDT is enabled
-        // DRVDT is enabled
+    for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) {
-        for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx)
+        if (oilVaporizationControl.drvdtActive(pvtRegionIdx)) {
             maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx) * timeStepSize;
+        }
     }
 }
 
@@ -177,16 +203,16 @@ updateMaxValues(const int episodeIdx,
                 const Scalar timeStepSize)
 {
     const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    if (this->drsdtActive(episodeIdx)) {
+    // DRSDT is enabled
-        // DRSDT is enabled
+    for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) {
-        for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) {
+        if (this->drsdtActive(episodeIdx, pvtRegionIdx)) {
             maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx) * timeStepSize;
         }
     }
 
-    if (this->drvdtActive(episodeIdx)) {
+       // DRVDT is enabled
-        // DRVDT is enabled
+    for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) {
-        for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) {
+        if (this->drvdtActive(episodeIdx, pvtRegionIdx)) {
             maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx) * timeStepSize;
         }
     }
@@ -207,6 +233,10 @@ drsdtcon(const unsigned elemIdx,
     // Output drsdt value for index 0
     episodeIdx = std::max(episodeIdx, 0);
     const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+
+    if (!oilVaporizationControl.drsdtConvective(pvtRegionIdx))
+        return 0;
+
     return oilVaporizationControl.getMaxDRSDT(pvtRegionIdx) * convectiveDrs_[elemIdx];
 }
 
@@ -218,12 +248,12 @@ maxGasDissolutionFactor(const unsigned timeIdx,
                         const int episodeIdx,
                         const int pvtRegionIdx) const
 {
-    if (!this->drsdtActive(episodeIdx) || maxDRs_[pvtRegionIdx] < 0.0) {
+    if (!this->drsdtActive(episodeIdx, pvtRegionIdx)) {
         return std::numeric_limits<Scalar>::max() / 2.0;
     }
 
     Scalar scaling = 1.0;
-    if (this->drsdtConvective(episodeIdx)) {
+    if (this->drsdtConvective(episodeIdx, pvtRegionIdx)) {
        scaling = convectiveDrs_[globalDofIdx];
     }
 
@@ -244,7 +274,7 @@ maxOilVaporizationFactor(const unsigned timeIdx,
                          const int episodeIdx,
                          const int pvtRegionIdx) const
 {
-    if (!this->drvdtActive(episodeIdx) || maxDRv_[pvtRegionIdx] < 0.0) {
+    if (!this->drvdtActive(episodeIdx, pvtRegionIdx)) {
         return std::numeric_limits<Scalar>::max() / 2.0;
     }
 

opm/simulators/flow/MixingRateControls.hpp

@@ -60,7 +60,10 @@ public:
     bool drsdtActive(int episodeIdx) const;
     bool drvdtActive(int episodeIdx) const;
     bool drsdtConvective(int episodeIdx) const;
-
+    
+    bool drsdtActive(int episodeIdx, std::size_t pvtRegionIdx) const;
+    bool drvdtActive(int episodeIdx, std::size_t pvtRegionIdx) const;
+    bool drsdtConvective(int episodeIdx, std::size_t pvtRegionIdx) const;
     /*!
      * \brief Returns the dynamic drsdt convective mixing value
      */
@@ -114,11 +117,11 @@ public:
                 const Scalar gravity,
                 const Scalar permZ,
                 const Scalar distZ,
-                const int pvtRegionIdx,
+                const int pvtRegionIdx)
-                const std::array<bool,3>& active)
     {
         const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-        if (active[0]) {
+
+        if (oilVaporizationControl.drsdtConvective(pvtRegionIdx)) {
             // This implements the convective DRSDT as described in
             // Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow
             // simulator" Submitted to TCCS 11, 2021
@@ -155,7 +158,7 @@ public:
                                          fs.pvtRegionIndex());
         }
 
-        if (active[1]) {
+        if (oilVaporizationControl.drsdtActive(pvtRegionIdx)) {
             const auto& fs = iq.fluidState();
 
             using FluidState = typename std::decay<decltype(fs)>::type;
@@ -171,7 +174,7 @@ public:
                 lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::infinity();
         }
 
-        if (active[2]) {
+        if (oilVaporizationControl.drvdtActive(pvtRegionIdx)) {
             const auto& fs = iq.fluidState();
             using FluidState = typename std::decay<decltype(fs)>::type;
             lastRv_[compressedDofIdx]