From b435c1b1bf25c68d63783b316d2c9e36beb0d2b8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Fri, 29 Nov 2024 13:39:28 +0100
Subject: [PATCH] For convective mixing, store the saturated dissolution
 factor.

---
 ...BlackOilIntensiveQuantitiesGlobalIndex.hpp |  7 +++
 .../blackoilconvectivemixingmodule.hh         | 58 ++++++++++++++++---
 .../blackoil/blackoilintensivequantities.hh   |  7 +++
 3 files changed, 63 insertions(+), 9 deletions(-)

diff --git a/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp b/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp
index 20f5b0d308e..e964b902a8a 100644
--- a/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp
+++ b/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp
@@ -49,6 +49,7 @@
 #include <opm/models/blackoil/blackoilenergymodules.hh>
 #include <opm/models/blackoil/blackoildiffusionmodule.hh>
 #include <opm/models/blackoil/blackoilmicpmodules.hh>
+#include <opm/models/blackoil/blackoilconvectivemixingmodule.hh>
 #include <opm/models/common/directionalmobility.hh>
 
 #include <opm/utility/CopyablePtr.hpp>
@@ -78,6 +79,7 @@ class BlackOilIntensiveQuantitiesGlobalIndex
     , public BlackOilBrineIntensiveQuantities<TypeTag>
     , public BlackOilEnergyIntensiveQuantitiesGlobalIndex<TypeTag>
     , public BlackOilMICPIntensiveQuantities<TypeTag>
+    , public BlackOilConvectiveMixingIntensiveQuantities<TypeTag>
 {
     using ParentType = GetPropType<TypeTag, Properties::DiscIntensiveQuantities>;
     using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
@@ -103,6 +105,7 @@ class BlackOilIntensiveQuantitiesGlobalIndex
     enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
+    enum { enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>() };
     enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
@@ -382,6 +385,10 @@ class BlackOilIntensiveQuantitiesGlobalIndex
 
         rockCompTransMultiplier_ = problem.template rockCompTransMultiplier<Evaluation>(*this, globalSpaceIdx);
 
+        if constexpr (enableConvectiveMixing) {
+            asImp_().updateSaturatedDissolutionFactor_();
+        }
+
 #ifndef NDEBUG
         // some safety checks in debug mode
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
diff --git a/opm/models/blackoil/blackoilconvectivemixingmodule.hh b/opm/models/blackoil/blackoilconvectivemixingmodule.hh
index ac46d3bfd64..1deac71b07e 100644
--- a/opm/models/blackoil/blackoilconvectivemixingmodule.hh
+++ b/opm/models/blackoil/blackoilconvectivemixingmodule.hh
@@ -289,15 +289,11 @@ public:
         const auto& liquidPhaseIdx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
             FluidSystem::waterPhaseIdx :
             FluidSystem::oilPhaseIdx;
-        const Evaluation SoMax = 0.0;
 
         //interiour
         const auto& t_in = intQuantsIn.fluidState().temperature(liquidPhaseIdx);
         const auto& p_in = intQuantsIn.fluidState().pressure(liquidPhaseIdx);
-        const auto& rssat_in = FluidSystem::saturatedDissolutionFactor(intQuantsIn.fluidState(),
-                                                            liquidPhaseIdx,
-                                                            intQuantsIn.pvtRegionIndex(),
-                                                            SoMax);
+        const auto& rssat_in = intQuantsIn.saturatedDissolutionFactor();
         const auto& salt_in = intQuantsIn.fluidState().saltSaturation();
 
         const auto bLiquidSatIn = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
@@ -315,10 +311,7 @@ public:
         //exteriour
         const auto t_ex = Opm::getValue(intQuantsEx.fluidState().temperature(liquidPhaseIdx));
         const auto p_ex = Opm::getValue(intQuantsEx.fluidState().pressure(liquidPhaseIdx));
-        const auto rssat_ex = Opm::getValue(FluidSystem::saturatedDissolutionFactor(intQuantsEx.fluidState(),
-                                                            liquidPhaseIdx,
-                                                            intQuantsEx.pvtRegionIndex(),
-                                                            SoMax));
+        const auto rssat_ex = Opm::getValue(intQuantsEx.saturatedDissolutionFactor());
         const auto salt_ex = Opm::getValue(intQuantsEx.fluidState().saltSaturation());
         const auto bLiquidSatEx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
             FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(), t_ex, p_ex, rssat_ex, salt_ex):
@@ -382,6 +375,53 @@ public:
     }
 };
 
+/*!
+ * \ingroup BlackOil
+ * \class Opm::BlackOilConvectiveMixingIntensiveQuantities
+ *
+ * \brief Provides the volumetric quantities required for the equations needed by the
+ *        convective mixing (DRSDTCON) model.
+ */
+template <class TypeTag, bool enableConvectiveMixingV = getPropValue<TypeTag, Properties::EnableConvectiveMixing>()>
+class BlackOilConvectiveMixingIntensiveQuantities
+{
+    using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+public:
+    /*!
+     * \brief Compute the intensive quantities needed to handle convective dissolution
+     *
+     */
+    void updateSaturatedDissolutionFactor_()
+    {
+        auto& fs = asImp_().fluidState();
+        const auto& liquidPhaseIdx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
+            FluidSystem::waterPhaseIdx :
+            FluidSystem::oilPhaseIdx;
+        const Evaluation SoMax = 0.0;
+        saturatedDissolutionFactor_ = FluidSystem::saturatedDissolutionFactor(fs,
+                                                                              liquidPhaseIdx,
+                                                                              asImp_().pvtRegionIndex(),
+                                                                              SoMax);
+    }
+
+    const Evaluation& saturatedDissolutionFactor() const
+    { return saturatedDissolutionFactor_; }
+
+protected:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    Evaluation saturatedDissolutionFactor_;
+};
+
+template <class TypeTag>
+class BlackOilConvectiveMixingIntensiveQuantities<TypeTag, false>
+{
+};
+
 }
 
 #endif
diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh
index 2de7bea2893..2d57b3e2b0f 100644
--- a/opm/models/blackoil/blackoilintensivequantities.hh
+++ b/opm/models/blackoil/blackoilintensivequantities.hh
@@ -38,6 +38,7 @@
 #include "blackoildiffusionmodule.hh"
 #include "blackoildispersionmodule.hh"
 #include "blackoilmicpmodules.hh"
+#include "blackoilconvectivemixingmodule.hh"
 
 #include <opm/common/TimingMacros.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
@@ -79,6 +80,7 @@ class BlackOilIntensiveQuantities
     , public BlackOilBrineIntensiveQuantities<TypeTag>
     , public BlackOilEnergyIntensiveQuantities<TypeTag>
     , public BlackOilMICPIntensiveQuantities<TypeTag>
+    , public BlackOilConvectiveMixingIntensiveQuantities<TypeTag>
 {
     using ParentType = GetPropType<TypeTag, Properties::DiscIntensiveQuantities>;
     using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
@@ -106,6 +108,7 @@ class BlackOilIntensiveQuantities
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
     enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
+    enum { enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>() };
     enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
@@ -497,6 +500,10 @@ public:
         // update the dispersion specific quantities of the intensive quantities
         DispersionIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
 
+        if constexpr (enableConvectiveMixing) {
+            asImp_().updateSaturatedDissolutionFactor_();
+        }
+
 #ifndef NDEBUG
         // some safety checks in debug mode
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {