From 7f37b16f3c66bea384c329a7890ceab4adda8a0e Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tosa@norceresearch.no>
Date: Wed, 26 Jun 2024 10:17:49 +0200
Subject: [PATCH] Add option to not use saturated tables

For co2storage and h2storage we dont have a concept of tables and should not spend time on
checking if we are at the saturated front
---
 .../fluidsystems/BlackOilFluidSystem.cpp      |  8 ++++
 .../fluidsystems/BlackOilFluidSystem.hpp      | 42 +++++++++++++------
 2 files changed, 38 insertions(+), 12 deletions(-)

diff --git a/opm/material/fluidsystems/BlackOilFluidSystem.cpp b/opm/material/fluidsystems/BlackOilFluidSystem.cpp
index af176d09dad..3d59ba6a753 100644
--- a/opm/material/fluidsystems/BlackOilFluidSystem.cpp
+++ b/opm/material/fluidsystems/BlackOilFluidSystem.cpp
@@ -159,6 +159,10 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
         }
     }
 
+    // For co2storage and h2storage we dont have a concept of tables and should not spend time on
+    // checking if we are at the saturated front
+    setUseSaturatedTables(!(eclState.runspec().h2Storage() || eclState.runspec().co2Storage()));
+
     setEnableDiffusion(eclState.getSimulationConfig().isDiffusive());
     if (enableDiffusion()) {
         const auto& diffCoeffTables = eclState.getTableManager().getDiffusionCoefficientTable();
@@ -220,6 +224,7 @@ void BlackOilFluidSystem<Scalar,IndexTraits>::
 initBegin(std::size_t numPvtRegions)
 {
     isInitialized_ = false;
+    useSaturatedTables_ = true;
 
     enableDissolvedGas_ = true;
     enableDissolvedGasInWater_ = false;
@@ -450,6 +455,9 @@ template <> std::vector<std::array<float, 9>> BlackOilFluidSystem<float, BlackOi
 template <> bool BlackOilFluidSystem<double, BlackOilDefaultIndexTraits>::isInitialized_ = false;
 template <> bool BlackOilFluidSystem<float, BlackOilDefaultIndexTraits>::isInitialized_ = false;
 
+template <> bool BlackOilFluidSystem<double, BlackOilDefaultIndexTraits>::useSaturatedTables_ = false;
+template <> bool BlackOilFluidSystem<float, BlackOilDefaultIndexTraits>::useSaturatedTables_ = false;
+
 // IMPORTANT: The following two lines must come after the template specializations above
 //    or else the static variable above will appear as undefined in the generated object file.
 template class BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>;
diff --git a/opm/material/fluidsystems/BlackOilFluidSystem.hpp b/opm/material/fluidsystems/BlackOilFluidSystem.hpp
index 5ce9b0e18b3..a7f7e158cbb 100644
--- a/opm/material/fluidsystems/BlackOilFluidSystem.hpp
+++ b/opm/material/fluidsystems/BlackOilFluidSystem.hpp
@@ -288,6 +288,13 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
     static void setEnableDiffusion(bool yesno)
     { enableDiffusion_ = yesno; }
 
+    /*!
+     * \brief Specify whether the saturated tables should be used
+     *
+     * By default, saturated tables are used
+     */
+    static void setUseSaturatedTables(bool yesno)
+    { useSaturatedTables_ = yesno; }
 
     /*!
      * \brief Set the pressure-volume-saturation (PVT) relations for the gas phase.
@@ -486,6 +493,14 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
     static bool enableDiffusion()
     { return enableDiffusion_; }
 
+    /*!
+     * \brief Returns whether the saturated tables should be used
+     *
+     * By default, saturated tables are used. If false the unsaturated tables are extrapolated
+     */
+    static bool useSaturatedTables()
+    { return useSaturatedTables_; }
+
     /*!
      * \brief Returns the density of a fluid phase at surface pressure [kg/m^3]
      *
@@ -751,7 +766,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
         case oilPhaseIdx: {
             if (enableDissolvedGas()) {
                 const auto& Rs = BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
-                if (fluidState.saturation(gasPhaseIdx) > 0.0
+                if (useSaturatedTables() && fluidState.saturation(gasPhaseIdx) > 0.0
                     && Rs >= (1.0 - 1e-10)*oilPvt_->saturatedGasDissolutionFactor(regionIdx, scalarValue(T), scalarValue(p)))
                 {
                     return oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
@@ -767,20 +782,20 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
             if (enableVaporizedOil() && enableVaporizedWater()) {
                  const auto& Rvw = BlackOil::template getRvw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                  const auto& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
-                 if (fluidState.saturation(waterPhaseIdx) > 0.0
+                 if (useSaturatedTables() && fluidState.saturation(waterPhaseIdx) > 0.0
                     && Rvw >= (1.0 - 1e-10)*gasPvt_->saturatedWaterVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p))
                     && fluidState.saturation(oilPhaseIdx) > 0.0
                     && Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                  { 
                     return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
                  } else {
-                     return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
+                    return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
                  }
             }
 
             if (enableVaporizedOil()) {
                 const auto& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
-                if (fluidState.saturation(oilPhaseIdx) > 0.0
+                if (useSaturatedTables() && fluidState.saturation(oilPhaseIdx) > 0.0
                     && Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                 {
                     return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
@@ -792,7 +807,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
 
             if (enableVaporizedWater()) { 
                 const auto& Rvw = BlackOil::template getRvw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
-                if (fluidState.saturation(waterPhaseIdx) > 0.0
+                if (useSaturatedTables() && fluidState.saturation(waterPhaseIdx) > 0.0
                     && Rvw >= (1.0 - 1e-10)*gasPvt_->saturatedWaterVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                 {
                     return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
@@ -811,7 +826,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
             const auto& saltConcentration = BlackOil::template getSaltConcentration_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
             if (enableDissolvedGasInWater()) {
                 const auto& Rsw = BlackOil::template getRsw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
-                if (fluidState.saturation(gasPhaseIdx) > 0.0
+                if (useSaturatedTables() && fluidState.saturation(gasPhaseIdx) > 0.0
                     && Rsw >= (1.0 - 1e-10)*waterPvt_->saturatedGasDissolutionFactor(regionIdx, scalarValue(T), scalarValue(p), scalarValue(saltConcentration)))
                 {
                     return waterPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p, saltConcentration);
@@ -995,7 +1010,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
         case oilPhaseIdx: {
             if (enableDissolvedGas()) {
                 const auto& Rs = BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
-                if (fluidState.saturation(gasPhaseIdx) > 0.0
+                if (useSaturatedTables() && fluidState.saturation(gasPhaseIdx) > 0.0
                     && Rs >= (1.0 - 1e-10)*oilPvt_->saturatedGasDissolutionFactor(regionIdx, scalarValue(T), scalarValue(p)))
                 {
                     return oilPvt_->saturatedViscosity(regionIdx, T, p);
@@ -1012,7 +1027,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
              if (enableVaporizedOil() && enableVaporizedWater()) {
                  const auto& Rvw = BlackOil::template getRvw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                  const auto& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
-                 if (fluidState.saturation(waterPhaseIdx) > 0.0
+                 if (useSaturatedTables() && fluidState.saturation(waterPhaseIdx) > 0.0
                     && Rvw >= (1.0 - 1e-10)*gasPvt_->saturatedWaterVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p))
                     && fluidState.saturation(oilPhaseIdx) > 0.0
                     && Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
@@ -1024,7 +1039,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
             }
             if (enableVaporizedOil()) {
                 const auto& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
-                if (fluidState.saturation(oilPhaseIdx) > 0.0
+                if (useSaturatedTables() && fluidState.saturation(oilPhaseIdx) > 0.0
                     && Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                 {
                     return gasPvt_->saturatedViscosity(regionIdx, T, p);
@@ -1035,7 +1050,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
             }
             if (enableVaporizedWater()) {
                 const auto& Rvw = BlackOil::template getRvw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
-                if (fluidState.saturation(waterPhaseIdx) > 0.0
+                if (useSaturatedTables() && fluidState.saturation(waterPhaseIdx) > 0.0
                     && Rvw >= (1.0 - 1e-10)*gasPvt_->saturatedWaterVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                 {
                     return gasPvt_->saturatedViscosity(regionIdx, T, p); 
@@ -1055,7 +1070,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
             const LhsEval& saltConcentration = BlackOil::template getSaltConcentration_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
             if (enableDissolvedGasInWater()) {
                 const auto& Rsw = BlackOil::template getRsw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
-                if (fluidState.saturation(gasPhaseIdx) > 0.0
+                if (useSaturatedTables() && fluidState.saturation(gasPhaseIdx) > 0.0
                     && Rsw >= (1.0 - 1e-10)*waterPvt_->saturatedGasDissolutionFactor(regionIdx, scalarValue(T), scalarValue(p), scalarValue(saltConcentration)))
                 {
                     return waterPvt_->saturatedViscosity(regionIdx, T, p, saltConcentration);
@@ -1541,6 +1556,8 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
     static std::array<short, numPhases> canonicalToActivePhaseIdx_;
 
     static bool isInitialized_;
+
+    static bool useSaturatedTables_;
 };
 
 template <typename T> using BOFS = BlackOilFluidSystem<T, BlackOilDefaultIndexTraits>;
@@ -1564,7 +1581,8 @@ template<> std::shared_ptr<WaterPvtMultiplexer<T>> BOFS<T>::waterPvt_; \
 template<> std::vector<std::array<T, 3>> BOFS<T>::referenceDensity_; \
 template<> std::vector<std::array<T, 3>> BOFS<T>::molarMass_; \
 template<> std::vector<std::array<T, 9>> BOFS<T>::diffusionCoefficients_; \
-template<> bool BOFS<T>::isInitialized_;
+template<> bool BOFS<T>::isInitialized_; \
+template<> bool BOFS<T>::useSaturatedTables_;
 
 DECLARE_INSTANCE(float)
 DECLARE_INSTANCE(double)