From a8aa12f0a2ec7e741fd6c07533ab3e8b82ed84d6 Mon Sep 17 00:00:00 2001
From: Paul Egberts <paul.egberts@tno.nl>
Date: Mon, 12 Sep 2022 16:02:35 +0200
Subject: [PATCH] added evaporated water dependency of gas density in equil
 calculations

---
 ebos/equil/initstateequil.cc | 66 ++++++++++++++++++++++++++++--------
 ebos/equil/initstateequil.hh |  6 ++--
 2 files changed, 55 insertions(+), 17 deletions(-)

diff --git a/ebos/equil/initstateequil.cc b/ebos/equil/initstateequil.cc
index 36e518aa2..35cb7ade7 100644
--- a/ebos/equil/initstateequil.cc
+++ b/ebos/equil/initstateequil.cc
@@ -325,30 +325,31 @@ density(const double depth,
     return rho;
 }
 
-// TODO: add RVW
-template<class FluidSystem, class RV>
-Gas<FluidSystem,RV>::
+template<class FluidSystem, class RV, class RVW>
+Gas<FluidSystem,RV,RVW>::
 Gas(const double temp,
         const RV& rv,
+        const RVW& rvw,
         const int pvtRegionIdx,
         const double normGrav)
         : temp_(temp)
         , rv_(rv)
+        , rvw_(rvw)
         , pvtRegionIdx_(pvtRegionIdx)
         , g_(normGrav)
 {
 }
 
-template<class FluidSystem, class RV>
-double Gas<FluidSystem,RV>::
+template<class FluidSystem, class RV, class RVW>
+double Gas<FluidSystem,RV,RVW>::
 operator()(const double depth,
            const double press) const
 {
     return this->density(depth, press) * g_;
 }
 
-template<class FluidSystem, class RV>
-double Gas<FluidSystem,RV>::
+template<class FluidSystem, class RV, class RVW>
+double Gas<FluidSystem,RV,RVW>::
 density(const double depth,
         const double press) const
 {
@@ -356,18 +357,53 @@ density(const double depth,
     if (FluidSystem::enableVaporizedOil())
         rv = rv_(depth, press, temp_);
 
+    double rvw = 0.0;
+    if (FluidSystem::enableVaporizedWater())
+        rvw = rvw_(depth, press, temp_);
+
     double bGas = 0.0;
-    if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)) {
-        bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
+
+    if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) {
+        if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)
+            && rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp_, press))
+        {
+            bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
+        } else {
+            bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, rvw);
+        }
+        double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
+        rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_)
+               + rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
+        return rho;
     }
-    else {
-        bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, 0.0/*=Rvw*/);
-    }
-    double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
-    if (FluidSystem::enableVaporizedOil()) {
+
+    if (FluidSystem::enableVaporizedOil()){
+        if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)) {
+            bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
+        } else {
+            bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, 0.0/*=rvw*/);
+        }
+        double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
         rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
+        return rho;
+     }
+
+    if (FluidSystem::enableVaporizedWater()){
+        if (rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp_, press)) {
+            bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
+        }
+        else {
+            bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0/*=rv*/, rvw);
+        }
+        double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
+        rho += rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
+        return rho;
     }
 
+    // immiscible gas
+    bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0/*=rv*/, 0.0/*=rvw*/);
+    double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
+
     return rho;
 }
 
@@ -1126,7 +1162,7 @@ makeGasPressure(const typename GPress::InitCond& ic,
                 const VSpan&                     span)
 {
     const auto drho = GasPressODE {
-        this->temperature_, reg.evaporationCalculator(),
+        this->temperature_, reg.evaporationCalculator(), reg.waterEvaporationCalculator(),
         reg.pvtIdx(), this->gravity_
     };
 
diff --git a/ebos/equil/initstateequil.hh b/ebos/equil/initstateequil.hh
index 8363a9118..30bca343e 100644
--- a/ebos/equil/initstateequil.hh
+++ b/ebos/equil/initstateequil.hh
@@ -129,12 +129,13 @@ private:
 };
 
 //TODO:add argument rvw
-template <class FluidSystem, class RV>
+template <class FluidSystem, class RV, class RVW>
 class Gas
 {
 public:
     Gas(const double temp,
         const RV& rv,
+        const RVW& rvw,
         const int pvtRegionIdx,
         const double normGrav);
 
@@ -144,6 +145,7 @@ public:
 private:
     const double temp_;
     const RV& rv_;
+    const RVW& rvw_;
     const int pvtRegionIdx_;
     const double g_;
 
@@ -277,7 +279,7 @@ private:
     >;
 
     using GasPressODE = PhasePressODE::Gas<
-        FluidSystem, typename Region::CalcEvaporation
+        FluidSystem, typename Region::CalcEvaporation, typename Region::CalcWaterEvaporation
     >;
 
     using WatPressODE = PhasePressODE::Water<FluidSystem>;