EclProblem: put drsdtcon with friends in separate class

CMakeLists_files.cmake

@@ -33,6 +33,7 @@ list (APPEND MAIN_SOURCE_FILES
   ebos/eclgenericvanguard.cc
   ebos/eclgenericwriter.cc
   ebos/eclinterregflows.cc
+  ebos/eclmixingratecontrols.cc
   ebos/eclsolutioncontainers.cc
   ebos/ecltransmissibility.cc
   ebos/equil/equilibrationhelpers.cc
@@ -386,6 +387,7 @@ list (APPEND PUBLIC_HEADER_FILES
   ebos/eclgenericwriter.hh
   ebos/eclgenericwriter_impl.hh
   ebos/eclinterregflows.hh
+  ebos/eclmixingratecontrols.hh
   ebos/eclmpiserializer.hh
   ebos/eclnewtonmethod.hh
   ebos/ecloutputblackoilmodule.hh

ebos/eclgenericproblem.hh

@@ -28,6 +28,7 @@
 #ifndef EWOMS_GENERIC_ECL_PROBLEM_HH
 #define EWOMS_GENERIC_ECL_PROBLEM_HH
 
+#include <ebos/eclmixingratecontrols.hh>
 #include <ebos/eclsolutioncontainers.hh>
 
 #include <opm/material/common/UniformXTabulated2DFunction.hpp>
@@ -275,8 +276,8 @@ public:
     { return maxFails_; }
 
     bool vapparsActive(int episodeIdx) const;
-    
+
-    int numPressurePointsEquil() const 
+    int numPressurePointsEquil() const
     { return numPressurePointsEquil_; }
 
     bool operator==(const EclGenericProblem& rhs) const;
@@ -291,22 +292,11 @@ public:
         serializer(overburdenPressure_);
         serializer(solventSaturation_);
         serializer(micp_);
-        serializer(lastRv_);
+        serializer(mixControls_);
-        serializer(maxDRv_);
-        serializer(convectiveDrs_);
-        serializer(lastRs_);
-        serializer(maxDRs_);
-        serializer(dRsDtOnlyFreeGas_);
     }
 
 protected:
-    bool drsdtActive_(int episodeIdx) const;
-    bool drvdtActive_(int episodeIdx) const;
-    bool drsdtConvective_(int episodeIdx) const;
-
     void initFluidSystem_();
-    void initDRSDT_(size_t numDof,
-                    int episodeIdx);
 
     /*!
      * \brief Always returns true. The ecl output writer takes care of the rest
@@ -383,13 +373,7 @@ protected:
     std::vector<Scalar> solventSaturation_;
     MICPSolutionContainer<Scalar> micp_;
 
-    std::vector<Scalar> lastRv_;
+    EclMixingRateControls<FluidSystem, Scalar> mixControls_;
-    std::vector<Scalar> maxDRv_;
-
-    std::vector<Scalar> convectiveDrs_;
-    std::vector<Scalar> lastRs_;
-    std::vector<Scalar> maxDRs_;
-    std::vector<bool> dRsDtOnlyFreeGas_; // apply the DRSDT rate limit only to cells that exhibit free gas
 
     // time stepping parameters
     bool enableTuning_;

ebos/eclgenericproblem_impl.hh

@@ -88,6 +88,7 @@ EclGenericProblem(const EclipseState& eclState,
     : eclState_(eclState)
     , schedule_(schedule)
     , gridView_(gridView)
+    , mixControls_(schedule)
 {
 }
 
@@ -106,12 +107,7 @@ serializationTestObject(const EclipseState& eclState,
     result.solventSaturation_ = {15.0};
     result.polymer_ = PolymerSolutionContainer<Scalar>::serializationTestObject();
     result.micp_ = MICPSolutionContainer<Scalar>::serializationTestObject();
-    result.lastRv_ = {21.0};
+    result.mixControls_ = EclMixingRateControls<FluidSystem,Scalar>::serializationTestObject(schedule);
-    result.maxDRv_ = {22.0, 23.0};
-    result.convectiveDrs_ = {24.0, 25.0, 26.0};
-    result.lastRs_ = {27.0};
-    result.maxDRs_ = {28.0};
-    result.dRsDtOnlyFreeGas_ = {false, true};
 
     return result;
 }
@@ -444,36 +440,6 @@ vapparsActive(int episodeIdx) const
     return (oilVaporizationControl.getType() == OilVaporizationProperties::OilVaporization::VAPPARS);
 }
 
-template<class GridView, class FluidSystem, class Scalar>
-bool EclGenericProblem<GridView,FluidSystem,Scalar>::
-drsdtActive_(int episodeIdx) const
-{
-    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
-                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
-    return (oilVaporizationControl.drsdtActive() && bothOilGasActive);
-}
-
-template<class GridView, class FluidSystem, class Scalar>
-bool EclGenericProblem<GridView,FluidSystem,Scalar>::
-drvdtActive_(int episodeIdx) const
-{
-    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
-                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
-    return (oilVaporizationControl.drvdtActive() && bothOilGasActive);
-}
-
-template<class GridView, class FluidSystem, class Scalar>
-bool EclGenericProblem<GridView,FluidSystem,Scalar>::
-drsdtConvective_(int episodeIdx) const
-{
-    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
-                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
-    return (oilVaporizationControl.drsdtConvective() && bothOilGasActive);
-}
-
 template<class GridView, class FluidSystem, class Scalar>
 bool EclGenericProblem<GridView,FluidSystem,Scalar>::
 beginEpisode_(bool enableExperiments,
@@ -539,16 +505,7 @@ beginTimeStep_(bool enableExperiments,
     }
 
     // update explicit quantities between timesteps.
-    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+    this->mixControls_.updateExplicitQuantities(episodeIdx, timeStepSize);
-    if (drsdtActive_(episodeIdx))
-        // DRSDT is enabled
-        for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx)
-            maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*timeStepSize;
-
-    if (drvdtActive_(episodeIdx))
-        // DRVDT is enabled
-        for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx)
-            maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*timeStepSize;
 }
 
 template<class GridView, class FluidSystem, class Scalar>
@@ -615,7 +572,7 @@ readBlackoilExtentionsInitialConditions_(size_t numDof,
         } else {
             micp_.calciteConcentration.resize(numDof, 0.0);
         }
-}
+    }
 }
 
 
@@ -663,15 +620,8 @@ template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 drsdtcon(unsigned elemIdx, int episodeIdx) const
 {
-    if (convectiveDrs_.empty())
+    return this->mixControls_.drsdtcon(elemIdx, episodeIdx,
-        return 0;
+                                       this->pvtRegionIndex(elemIdx));
-
-    // The episode index is set to -1 in the initialization phase.
-    // Output drsdt value for index 0
-    episodeIdx = std::max(episodeIdx, 0);
-    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    int pvtRegionIdx = pvtRegionIndex(elemIdx);
-    return oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*convectiveDrs_[elemIdx];
 }
 
 template<class GridView, class FluidSystem, class Scalar>
@@ -802,48 +752,17 @@ maxPolymerAdsorption(unsigned elemIdx) const
     return polymer_.maxAdsorption[elemIdx];
 }
 
-template<class GridView, class FluidSystem, class Scalar>
-void EclGenericProblem<GridView,FluidSystem,Scalar>::
-initDRSDT_(size_t numDof,
-           int episodeIdx)
-{
-    // deal with DRSDT
-    unsigned ntpvt = eclState_.runspec().tabdims().getNumPVTTables();
-    //TODO We may want to only allocate these properties only if active.
-    //But since they may be activated at later time we need some more
-    //intrastructure to handle it
-    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-        maxDRv_.resize(ntpvt, 1e30);
-        lastRv_.resize(numDof, 0.0);
-        maxDRs_.resize(ntpvt, 1e30);
-        dRsDtOnlyFreeGas_.resize(ntpvt, false);
-        lastRs_.resize(numDof, 0.0);
-        maxDRv_.resize(ntpvt, 1e30);
-        lastRv_.resize(numDof, 0.0);
-        maxOilSaturation_.resize(numDof, 0.0);
-        if (drsdtConvective_(episodeIdx)) {
-            convectiveDrs_.resize(numDof, 1.0);
-        }
-    }
-}
-
 template<class GridView, class FluidSystem, class Scalar>
 bool EclGenericProblem<GridView,FluidSystem,Scalar>::
 operator==(const EclGenericProblem& rhs) const
 {
-    return this->maxOilSaturation_ == rhs.maxOilSaturation_ &&
+    return this->maxWaterSaturation_ == rhs.maxWaterSaturation_ &&
-           this->maxWaterSaturation_ == rhs.maxWaterSaturation_ &&
            this->minOilPressure_ == rhs.minOilPressure_ &&
            this->overburdenPressure_ == rhs.overburdenPressure_ &&
            this->solventSaturation_ == rhs.solventSaturation_ &&
            this->polymer_ == rhs.polymer_ &&
            this->micp_ == rhs.micp_ &&
-           this->lastRv_ == rhs.lastRv_ &&
+           this->mixControls_ == rhs.mixControls_;
-           this->maxDRv_ == rhs.maxDRv_ &&
-           this->convectiveDrs_ == rhs.convectiveDrs_ &&
-           this->lastRs_ == rhs.lastRs_ &&
-           this->maxDRs_ == rhs.maxDRs_ &&
-           this->dRsDtOnlyFreeGas_ == rhs.dRsDtOnlyFreeGas_;
 }
 
 } // namespace Opm

ebos/eclmixingratecontrols.cc

@@ -0,0 +1,299 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::EclProblem
+ */
+
+#include <config.h>
+#include <ebos/eclmixingratecontrols.hh>
+
+#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
+
+#include <algorithm>
+#include <limits>
+
+namespace Opm {
+
+template<class FluidSystem, class Scalar>
+EclMixingRateControls<FluidSystem,Scalar>::
+EclMixingRateControls(const Schedule& schedule)
+    : schedule_(schedule)
+{}
+
+template<class FluidSystem, class Scalar>
+EclMixingRateControls<FluidSystem,Scalar>::
+EclMixingRateControls(const EclMixingRateControls& rhs)
+    : schedule_(rhs.schedule_)
+{
+    *this = rhs;
+}
+
+template<class FluidSystem, class Scalar>
+EclMixingRateControls<FluidSystem,Scalar>
+EclMixingRateControls<FluidSystem,Scalar>::
+serializationTestObject(const Schedule& schedule)
+{
+    EclMixingRateControls<FluidSystem,Scalar> result(schedule);
+    result.lastRv_ = {21.0};
+    result.maxDRv_ = {22.0, 23.0};
+    result.convectiveDrs_ = {24.0, 25.0, 26.0};
+    result.lastRs_ = {27.0};
+    result.maxDRs_ = {28.0};
+    result.dRsDtOnlyFreeGas_ = {false, true};
+
+    return result;
+}
+
+template<class FluidSystem, class Scalar>
+bool EclMixingRateControls<FluidSystem,Scalar>::
+operator==(const EclMixingRateControls& rhs) const
+{
+    return this->lastRv_ == rhs.lastRv_ &&
+           this->maxDRv_ == rhs.maxDRv_ &&
+           this->convectiveDrs_ == rhs.convectiveDrs_ &&
+           this->lastRs_ == rhs.lastRs_ &&
+           this->maxDRs_ == rhs.maxDRs_ &&
+           this->dRsDtOnlyFreeGas_ == rhs.dRsDtOnlyFreeGas_;
+}
+
+template<class FluidSystem, class Scalar>
+EclMixingRateControls<FluidSystem,Scalar>&
+EclMixingRateControls<FluidSystem,Scalar>::
+operator=(const EclMixingRateControls& rhs)
+{
+    this->lastRv_ = rhs.lastRv_;
+    this->maxDRv_ = rhs.maxDRv_;
+    this->convectiveDrs_ = rhs.convectiveDrs_;
+    this->lastRs_ = rhs.lastRs_;
+    this->maxDRs_ = rhs.maxDRs_;
+    this->dRsDtOnlyFreeGas_ = rhs.dRsDtOnlyFreeGas_;
+
+    return *this;
+}
+
+template<class FluidSystem, class Scalar>
+void EclMixingRateControls<FluidSystem,Scalar>::
+init(std::size_t numDof, int episodeIdx, const unsigned ntpvt)
+{
+    // deal with DRSDT
+    //TODO We may want to only allocate these properties only if active.
+    //But since they may be activated at later time we need some more
+    //intrastructure to handle it
+    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
+        FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+        maxDRv_.resize(ntpvt, 1e30);
+        lastRv_.resize(numDof, 0.0);
+        maxDRs_.resize(ntpvt, 1e30);
+        dRsDtOnlyFreeGas_.resize(ntpvt, false);
+        lastRs_.resize(numDof, 0.0);
+        maxDRv_.resize(ntpvt, 1e30);
+        lastRv_.resize(numDof, 0.0);
+        if (this->drsdtConvective(episodeIdx)) {
+            convectiveDrs_.resize(numDof, 1.0);
+        }
+    }
+}
+
+template<class FluidSystem, class Scalar>
+bool EclMixingRateControls<FluidSystem,Scalar>::
+drsdtActive(int episodeIdx) const
+{
+    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
+                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
+    return (oilVaporizationControl.drsdtActive() && bothOilGasActive);
+}
+
+template<class FluidSystem, class Scalar>
+bool EclMixingRateControls<FluidSystem,Scalar>::
+drvdtActive(int episodeIdx) const
+{
+    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
+                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
+    return (oilVaporizationControl.drvdtActive() && bothOilGasActive);
+}
+
+template<class FluidSystem, class Scalar>
+bool EclMixingRateControls<FluidSystem,Scalar>::
+drsdtConvective(int episodeIdx) const
+{
+    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
+                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
+    return (oilVaporizationControl.drsdtConvective() && bothOilGasActive);
+}
+
+template<class FluidSystem, class Scalar>
+void EclMixingRateControls<FluidSystem,Scalar>::
+updateExplicitQuantities(const int episodeIdx,
+                         const Scalar timeStepSize)
+{
+    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+    if (this->drsdtActive(episodeIdx)) {
+        // DRSDT is enabled
+        for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx)
+            maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx) * timeStepSize;
+    }
+
+    if (this->drvdtActive(episodeIdx)) {
+        // DRVDT is enabled
+        for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx)
+            maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx) * timeStepSize;
+    }
+}
+
+template<class FluidSystem, class Scalar>
+void EclMixingRateControls<FluidSystem,Scalar>::
+updateLastValues(const unsigned elemIdx,
+                 const Scalar Rs,
+                 const Scalar Rv)
+{
+    if (!lastRs_.empty()) {
+        lastRs_[elemIdx] = Rs;
+    }
+
+    if (!lastRv_.empty()) {
+        lastRv_[elemIdx] = Rv;
+    }
+}
+
+template<class FluidSystem, class Scalar>
+void EclMixingRateControls<FluidSystem,Scalar>::
+updateMaxValues(const int episodeIdx,
+                const Scalar timeStepSize)
+{
+    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+    if (this->drsdtActive(episodeIdx)) {
+        // DRSDT is enabled
+        for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) {
+            maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx) * timeStepSize;
+        }
+    }
+
+    if (this->drvdtActive(episodeIdx)) {
+        // DRVDT is enabled
+        for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) {
+            maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx) * timeStepSize;
+        }
+    }
+}
+
+template<class FluidSystem, class Scalar>
+Scalar EclMixingRateControls<FluidSystem,Scalar>::
+drsdtcon(const unsigned elemIdx,
+         int episodeIdx,
+         const int pvtRegionIdx) const
+{
+    if (convectiveDrs_.empty()) {
+        return 0;
+    }
+
+    // The episode index is set to -1 in the initialization phase.
+    // Output drsdt value for index 0
+    episodeIdx = std::max(episodeIdx, 0);
+    const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+    return oilVaporizationControl.getMaxDRSDT(pvtRegionIdx) * convectiveDrs_[elemIdx];
+}
+
+template<class FluidSystem, class Scalar>
+Scalar EclMixingRateControls<FluidSystem,Scalar>::
+maxGasDissolutionFactor(const unsigned timeIdx,
+                        const unsigned globalDofIdx,
+                        const int episodeIdx,
+                        const int pvtRegionIdx) const
+{
+    if (!this->drsdtActive(episodeIdx) || maxDRs_[pvtRegionIdx] < 0.0) {
+        return std::numeric_limits<Scalar>::max() / 2.0;
+    }
+
+    Scalar scaling = 1.0;
+    if (this->drsdtConvective(episodeIdx)) {
+       scaling = convectiveDrs_[globalDofIdx];
+    }
+
+    // this is a bit hacky because it assumes that a time discretization with only
+    // two time indices is used.
+    if (timeIdx == 0) {
+        return lastRs_[globalDofIdx] + maxDRs_[pvtRegionIdx] * scaling;
+    } else {
+        return lastRs_[globalDofIdx];
+    }
+}
+
+template<class FluidSystem, class Scalar>
+Scalar EclMixingRateControls<FluidSystem,Scalar>::
+maxOilVaporizationFactor(const unsigned timeIdx,
+                         const unsigned globalDofIdx,
+                         const int episodeIdx,
+                         const int pvtRegionIdx) const
+{
+    if (!this->drvdtActive(episodeIdx) || maxDRv_[pvtRegionIdx] < 0.0) {
+        return std::numeric_limits<Scalar>::max() / 2.0;
+    }
+
+    // this is a bit hacky because it assumes that a time discretization with only
+    // two time indices is used.
+    if (timeIdx == 0) {
+        return lastRv_[globalDofIdx] + maxDRv_[pvtRegionIdx];
+    } else {
+        return lastRv_[globalDofIdx];
+    }
+}
+template<class FluidSystem, class Scalar>
+void EclMixingRateControls<FluidSystem,Scalar>::
+updateConvectiveDRsDt_(const unsigned compressedDofIdx,
+                       const Scalar t,
+                       const Scalar p,
+                       const Scalar rs,
+                       const Scalar so,
+                       const Scalar poro,
+                       const Scalar permz,
+                       const Scalar distZ,
+                       const Scalar gravity,
+                       const int pvtRegionIndex)
+{
+    const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p);
+    const Scalar saturatedInvB
+        = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p);
+    const Scalar rsZero = 0.0;
+    const Scalar pureDensity
+        = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero)
+        * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex);
+    const Scalar saturatedDensity = saturatedInvB
+        * (FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex)
+           + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex));
+    const Scalar deltaDensity = saturatedDensity - pureDensity;
+    const Scalar visc = FluidSystem::oilPvt().viscosity(pvtRegionIndex, t, p, rs);
+    // Note that for so = 0 this gives no limits (inf) for the dissolution rate
+    // Also we restrict the effect of convective mixing to positive density differences
+    // i.e. we only allow for fingers moving downward
+    convectiveDrs_[compressedDofIdx]
+        = permz * rssat * max(0.0, deltaDensity) * gravity / (so * visc * distZ * poro);
+}
+
+template class EclMixingRateControls<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>, double>;
+
+} // namespace Opm

ebos/eclmixingratecontrols.hh

@@ -0,0 +1,183 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::EclProblem
+ */
+#ifndef ECL_MIXING_RATE_CONTROLS_HH
+#define ECL_MIXING_RATE_CONTROLS_HH
+
+#include <opm/input/eclipse/Schedule/Schedule.hpp>
+
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+#include <opm/material/common/MathToolbox.hpp>
+
+#include <limits>
+#include <vector>
+
+namespace Opm {
+
+class EclipseState;
+
+//! \brief Class handling mixing rate controls for an EclProblem.
+template<class FluidSystem, class Scalar>
+class EclMixingRateControls {
+public:
+    EclMixingRateControls(const Schedule& schedule);
+    EclMixingRateControls(const EclMixingRateControls& rhs);
+
+    static EclMixingRateControls serializationTestObject(const Schedule& schedule);
+
+    bool operator==(const EclMixingRateControls& rhs) const;
+    EclMixingRateControls& operator=(const EclMixingRateControls& rhs);
+
+    void init(std::size_t numDof, int episodeIdx, const unsigned ntpvt);
+
+    bool drsdtActive(int episodeIdx) const;
+    bool drvdtActive(int episodeIdx) const;
+    bool drsdtConvective(int episodeIdx) const;
+
+    /*!
+     * \brief Returns the dynamic drsdt convective mixing value
+     */
+    Scalar drsdtcon(const unsigned elemIdx,
+                    int episodeIdx,
+                    const int pvtRegionIdx) const;
+
+    /*!
+     * \brief Returns the maximum value of the gas dissolution factor at the current time
+     *        for a given degree of freedom.
+     */
+    Scalar maxGasDissolutionFactor(unsigned timeIdx,
+                                   unsigned globalDofIdx,
+                                   const int episodeIdx,
+                                   const int pvtRegionIdx) const;
+
+    /*!
+     * \brief Returns the maximum value of the oil vaporization factor at the current
+     *        time for a given degree of freedom.
+     */
+    Scalar maxOilVaporizationFactor(const unsigned timeIdx,
+                                    const unsigned globalDofIdx,
+                                    const int episodeIdx,
+                                    const int pvtRegionIdx) const;
+
+    void updateExplicitQuantities(const int episodeIdx,
+                                  const Scalar timeStepSize);
+
+    void updateLastValues(const unsigned elemIdx,
+                          const Scalar Rs,
+                          const Scalar Rv);
+
+    void updateMaxValues(const int episodeIdx,
+                         const Scalar timeStepSize);
+
+    template<class Serializer>
+    void serializeOp(Serializer& serializer)
+    {
+        serializer(lastRv_);
+        serializer(maxDRv_);
+        serializer(convectiveDrs_);
+        serializer(lastRs_);
+        serializer(maxDRs_);
+        serializer(dRsDtOnlyFreeGas_);
+    }
+
+    template<class IntensiveQuantities>
+    void update(unsigned compressedDofIdx,
+                const IntensiveQuantities& iq,
+                const int episodeIdx,
+                const Scalar gravity,
+                const Scalar permZ,
+                const Scalar distZ,
+                const int pvtRegionIdx,
+                const std::array<bool,3>& active)
+    {
+        if (active[0]) {
+            // 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
+            const auto& fs = iq.fluidState();
+            this->updateConvectiveDRsDt_(compressedDofIdx,
+                                         getValue(fs.temperature(FluidSystem::oilPhaseIdx)),
+                                         getValue(fs.pressure(FluidSystem::oilPhaseIdx)),
+                                         getValue(fs.Rs()),
+                                         getValue(fs.saturation(FluidSystem::oilPhaseIdx)),
+                                         getValue(iq.porosity()),
+                                         permZ,
+                                         distZ,
+                                         gravity,
+                                         fs.pvtRegionIndex());
+        }
+
+        if (active[1]) {
+            const auto& fs = iq.fluidState();
+
+            using FluidState = typename std::decay<decltype(fs)>::type;
+
+            const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
+            constexpr Scalar freeGasMinSaturation_ = 1e-7;
+            if (oilVaporizationControl.getOption(pvtRegionIdx) ||
+                fs.saturation(FluidSystem::gasPhaseIdx) > freeGasMinSaturation_) {
+                lastRs_[compressedDofIdx]
+                    = BlackOil::template getRs_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
+            }
+            else
+                lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::infinity();
+        }
+
+        if (active[2]) {
+            const auto& fs = iq.fluidState();
+            using FluidState = typename std::decay<decltype(fs)>::type;
+            lastRv_[compressedDofIdx]
+                = BlackOil::template getRv_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
+        }
+    }
+
+private:
+    void updateConvectiveDRsDt_(const unsigned compressedDofIdx,
+                                const Scalar t,
+                                const Scalar p,
+                                const Scalar rs,
+                                const Scalar so,
+                                const Scalar poro,
+                                const Scalar permz,
+                                const Scalar distZ,
+                                const Scalar gravity,
+                                const int pvtRegionIndex);
+
+    std::vector<Scalar> lastRv_;
+    std::vector<Scalar> maxDRv_;
+
+    std::vector<Scalar> convectiveDrs_;
+    std::vector<Scalar> lastRs_;
+    std::vector<Scalar> maxDRs_;
+    std::vector<bool> dRsDtOnlyFreeGas_; // apply the DRSDT rate limit only to cells that exhibit free gas
+
+    const Schedule& schedule_;
+};
+
+} // namespace Opm
+
+#endif

ebos/eclproblem.hh

@@ -375,7 +375,14 @@ public:
         this->initFluidSystem_();
 
         // deal with DRSDT
-        this->initDRSDT_(this->model().numGridDof(), this->episodeIndex());
+        this->mixControls_.init(this->model().numGridDof(),
+                                this->episodeIndex(),
+                                eclState.runspec().tabdims().getNumPVTTables());
+
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
+            FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
+        }
 
         this->readRockParameters_(simulator.vanguard().cellCenterDepths(),
                                   [&simulator](const unsigned idx)
@@ -1166,22 +1173,9 @@ public:
      */
     Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const
     {
-        int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
+        return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx,
-        int episodeIdx = this->episodeIndex();
+                                                          this->episodeIndex(),
-        if (!this->drsdtActive_(episodeIdx) || this->maxDRs_[pvtRegionIdx] < 0.0)
+                                                          this->pvtRegionIndex(globalDofIdx));
-            return std::numeric_limits<Scalar>::max()/2.0;
-
-        Scalar scaling = 1.0;
-        if (this->drsdtConvective_(episodeIdx)) {
-           scaling = this->convectiveDrs_[globalDofIdx];
-        }
-
-        // this is a bit hacky because it assumes that a time discretization with only
-        // two time indices is used.
-        if (timeIdx == 0)
-            return this->lastRs_[globalDofIdx] + this->maxDRs_[pvtRegionIdx] * scaling;
-        else
-            return this->lastRs_[globalDofIdx];
     }
 
     /*!
@@ -1190,17 +1184,9 @@ public:
      */
     Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const
     {
-        int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
+        return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx,
-        int episodeIdx = this->episodeIndex();
+                                                           this->episodeIndex(),
-        if (!this->drvdtActive_(episodeIdx) || this->maxDRv_[pvtRegionIdx] < 0.0)
+                                                           this->pvtRegionIndex(globalDofIdx));
-            return std::numeric_limits<Scalar>::max()/2.0;
-
-        // this is a bit hacky because it assumes that a time discretization with only
-        // two time indices is used.
-        if (timeIdx == 0)
-            return this->lastRv_[globalDofIdx] + this->maxDRv_[pvtRegionIdx];
-        else
-            return this->lastRv_[globalDofIdx];
     }
 
     /*!
@@ -1214,8 +1200,8 @@ public:
     bool recycleFirstIterationStorage() const
     {
         int episodeIdx = this->episodeIndex();
-        return !this->drsdtActive_(episodeIdx) &&
+        return !this->mixControls_.drsdtActive(episodeIdx) &&
-               !this->drvdtActive_(episodeIdx) &&
+               !this->mixControls_.drvdtActive(episodeIdx) &&
                this->rockCompPoroMultWc_.empty() &&
                this->rockCompPoroMult_.empty();
     }
@@ -1712,82 +1698,32 @@ protected:
         // update the "last Rs" values for all elements, including the ones in the ghost
         // and overlap regions
         int episodeIdx = this->episodeIndex();
-        std::array<bool,3> active{this->drsdtConvective_(episodeIdx),
+        std::array<bool,3> active{this->mixControls_.drsdtConvective(episodeIdx),
-                                  this->drsdtActive_(episodeIdx),
+                                  this->mixControls_.drsdtActive(episodeIdx),
-                                  this->drvdtActive_(episodeIdx)};
+                                  this->mixControls_.drvdtActive(episodeIdx)};
-        if (!active[0] && !active[1] && !active[2])
+        if (!active[0] && !active[1] && !active[2]) {
-          return;
+            return;
+        }
 
         this->updateProperty_("EclProblem::updateCompositionChangeLimits_()) failed:",
-                              [this,episodeIdx,active](unsigned compressedDofIdx, const IntensiveQuantities& iq)
+                              [this,episodeIdx,active](unsigned compressedDofIdx,
+                                                       const IntensiveQuantities& iq)
                               {
-                                  this->updateCompositionChangeLimits_(compressedDofIdx,
+                                  const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx);
-                                                                       iq,
+                                  const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0;
-                                                                       episodeIdx,
+                                  const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
-                                                                       active);
+                                  this->mixControls_.update(compressedDofIdx,
+                                                            iq,
+                                                            episodeIdx,
+                                                            this->gravity_[dim - 1],
+                                                            perm[dim - 1][dim - 1],
+                                                            distZ,
+                                                            pvtRegionIdx,
+                                                            active);
                               }
             );
     }
 
-    void updateCompositionChangeLimits_(unsigned compressedDofIdx, const IntensiveQuantities& iq,int episodeIdx, const std::array<bool,3>& active)
-    {
-        auto& simulator = this->simulator();
-        auto& vanguard = simulator.vanguard();
-        if (active[0]) {
-            // 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
-            const Scalar g = this->gravity_[dim - 1];
-            const DimMatrix& perm = intrinsicPermeability(compressedDofIdx);
-            const Scalar permz = perm[dim - 1][dim - 1]; // The Z permeability
-            const Scalar distZ = vanguard.cellThickness(compressedDofIdx);
-            const auto& fs = iq.fluidState();
-            const Scalar t = getValue(fs.temperature(FluidSystem::oilPhaseIdx));
-            const Scalar p = getValue(fs.pressure(FluidSystem::oilPhaseIdx));
-            const Scalar so = getValue(fs.saturation(FluidSystem::oilPhaseIdx));
-            const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), t, p);
-            const Scalar saturatedInvB
-                = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), t, p);
-            const Scalar rsZero = 0.0;
-            const Scalar pureDensity
-                = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), t, p, rsZero)
-                * FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex());
-            const Scalar saturatedDensity = saturatedInvB
-                * (FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex())
-                   + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, fs.pvtRegionIndex()));
-            const Scalar deltaDensity = saturatedDensity - pureDensity;
-            const Scalar rs = getValue(fs.Rs());
-            const Scalar visc = FluidSystem::oilPvt().viscosity(fs.pvtRegionIndex(), t, p, rs);
-            const Scalar poro = getValue(iq.porosity());
-            // Note that for so = 0 this gives no limits (inf) for the dissolution rate
-            // Also we restrict the effect of convective mixing to positive density differences
-            // i.e. we only allow for fingers moving downward
-            this->convectiveDrs_[compressedDofIdx]
-                = permz * rssat * max(0.0, deltaDensity) * g / (so * visc * distZ * poro);
-        }
-
-        if (active[1]) {
-            const auto& fs = iq.fluidState();
-
-            using FluidState = typename std::decay<decltype(fs)>::type;
-
-            int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
-            const auto& oilVaporizationControl = vanguard.schedule()[episodeIdx].oilvap();
-            if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_)
-                this->lastRs_[compressedDofIdx]
-                    = BlackOil::template getRs_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
-            else
-                this->lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::infinity();
-        }
-
-        if (active[2]) {
-            const auto& fs = iq.fluidState();
-            using FluidState = typename std::decay<decltype(fs)>::type;
-            this->lastRv_[compressedDofIdx]
-                = BlackOil::template getRv_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
-        }
-    }
-
     bool updateMaxOilSaturation_()
     {
         OPM_TIMEBLOCK(updateMaxOilSaturation);
@@ -2065,13 +2001,7 @@ protected:
                     this->solventSaturation_[elemIdx] = ssol;
             }
 
-            if (! this->lastRs_.empty()) {
+            this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
-                this->lastRs_[elemIdx] = elemFluidState.Rs();
-            }
-
-            if (! this->lastRv_.empty()) {
-                this->lastRv_[elemIdx] = elemFluidState.Rv();
-            }
 
             if constexpr (enablePolymer)
                  this->polymer_.concentration[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx);
@@ -2086,16 +2016,7 @@ protected:
         }
 
         const int episodeIdx = this->episodeIndex();
-        const auto& oilVaporizationControl = simulator.vanguard().schedule()[episodeIdx].oilvap();
+        this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize());
-        if (this->drsdtActive_(episodeIdx))
-            // DRSDT is enabled
-            for (size_t pvtRegionIdx = 0; pvtRegionIdx < this->maxDRs_.size(); ++pvtRegionIdx)
-                this->maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*simulator.timeStepSize();
-
-        if (this->drvdtActive_(episodeIdx))
-            // DRVDT is enabled
-            for (size_t pvtRegionIdx = 0; pvtRegionIdx < this->maxDRv_.size(); ++pvtRegionIdx)
-                this->maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*simulator.timeStepSize();
 
         // assign the restart solution to the current solution. note that we still need
         // to compute real initial solution after this because the initial fluid states
@@ -2554,8 +2475,6 @@ private:
 
     std::vector<InitialFluidState> initialFluidStates_;
 
-    constexpr static Scalar freeGasMinSaturation_ = 1e-7;
-
     bool enableDriftCompensation_;
     GlobalEqVector drift_;