Merge pull request #4292 from totto82/addRsw

Add support for gas in water

CMakeLists.txt

@@ -372,7 +372,8 @@ set(FLOW_MODELS blackoil brine energy extbo foam gasoil gaswater
                 oilwater_polymer_injectivity micp polymer solvent
                 gasoil_energy brine_saltprecipitation
                 gaswater_saltprec_vapwat brine_precsalt_vapwat
-                blackoil_legacyassembly gasoildiffuse)
+                blackoil_legacyassembly gasoildiffuse gaswater_dissolution
+                gaswater_dissolution_diffuse)
 set(FLOW_VARIANT_MODELS brine_energy onephase onephase_energy)
 
 set(FLOW_TGTS)

ebos/eclequilinitializer.hh

@@ -77,6 +77,8 @@ class EclEquilInitializer
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
     enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
     enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
+    enum { has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
+
 
 public:
     // NB: setting the enableEnergy argument to true enables storage of enthalpy and
@@ -89,6 +91,7 @@ public:
                                                 enableEvaporation,
                                                 enableBrine,
                                                 enableSaltPrecipitation,
+                                                has_disgas_in_water,
                                                 Indices::numPhases>;
 
 

ebos/eclgenericoutputblackoilmodule.cc

@@ -685,6 +685,7 @@ assignToSolution(data::Solution& sol)
         {"PCOG",     UnitSystem::measure::pressure,  data::TargetType::RESTART_SOLUTION,      pcog_},
         {"PRES_OVB", UnitSystem::measure::pressure,  data::TargetType::RESTART_SOLUTION,      overburdenPressure_},
         {"RS",       UnitSystem::measure::gas_oil_ratio, data::TargetType::RESTART_SOLUTION,  rs_},
+        {"RSW",      UnitSystem::measure::gas_oil_ratio, data::TargetType::RESTART_SOLUTION,  rsw_},
         {"RSSAT",    UnitSystem::measure::gas_oil_ratio, data::TargetType::RESTART_AUXILIARY, gasDissolutionFactor_},
         {"RV",       UnitSystem::measure::oil_gas_ratio, data::TargetType::RESTART_SOLUTION,  rv_},
         {"RVSAT",    UnitSystem::measure::oil_gas_ratio, data::TargetType::RESTART_AUXILIARY, oilVaporizationFactor_},
@@ -791,6 +792,8 @@ setRestart(const data::Solution& sol,
         temperature_[elemIdx] = sol.data("TEMP")[globalDofIndex];
     if (!rs_.empty() && sol.has("RS"))
         rs_[elemIdx] = sol.data("RS")[globalDofIndex];
+    if (!rsw_.empty() && sol.has("RSW"))
+        rsw_[elemIdx] = sol.data("RSW")[globalDofIndex];
     if (!rv_.empty() && sol.has("RV"))
         rv_[elemIdx] = sol.data("RV")[globalDofIndex];
     if (!rvw_.empty() && sol.has("RVW"))
@@ -985,6 +988,10 @@ doAllocBuffers(unsigned bufferSize,
         rs_.resize(bufferSize, 0.0);
         rstKeywords["RS"] = 0;
     }
+    if (FluidSystem::enableDissolvedGasInWater()) {
+        rsw_.resize(bufferSize, 0.0);
+        rstKeywords["RSW"] = 0;
+    }
     if (FluidSystem::enableVaporizedOil()) {
         rv_.resize(bufferSize, 0.0);
         rstKeywords["RV"] = 0;

ebos/eclgenericoutputblackoilmodule.hh

@@ -414,6 +414,7 @@ protected:
     ScalarBuffer oilPressure_;
     ScalarBuffer temperature_;
     ScalarBuffer rs_;
+    ScalarBuffer rsw_;
     ScalarBuffer rv_;
     ScalarBuffer rvw_;
     ScalarBuffer overburdenPressure_;

ebos/ecloutputblackoilmodule.hh

@@ -285,6 +285,10 @@ public:
                 this->rs_[globalDofIdx] = getValue(fs.Rs());
                 Valgrind::CheckDefined(this->rs_[globalDofIdx]);
             }
+            if (!this->rsw_.empty()) {
+                this->rsw_[globalDofIdx] = getValue(fs.Rsw());
+                Valgrind::CheckDefined(this->rsw_[globalDofIdx]);
+            }
 
             if (!this->rv_.empty()) {
                 this->rv_[globalDofIdx] = getValue(fs.Rv());
@@ -498,6 +502,9 @@ public:
                 if (!this->rs_.empty())
                     this->rs_[globalDofIdx] = fsInitial.Rs();
 
+                if (!this->rsw_.empty())
+                    this->rsw_[globalDofIdx] = fsInitial.Rsw();
+
                 if (!this->rvw_.empty())
                     this->rvw_[globalDofIdx] = fsInitial.Rvw();
 
@@ -823,6 +830,8 @@ public:
             fs.setTemperature(this->temperature_[elemIdx]);
         if (!this->rs_.empty())
            fs.setRs(this->rs_[elemIdx]);
+        if (!this->rsw_.empty())
+           fs.setRsw(this->rsw_[elemIdx]);
         if (!this->rv_.empty())
            fs.setRv(this->rv_[elemIdx]);
         if (!this->rvw_.empty())

ebos/equil/initstateequil.cc

@@ -280,7 +280,7 @@ density(const double depth,
 {
     // The initializing algorithm can give depths outside the range due to numerical noise i.e. we extrapolate
     double saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true);
-    double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, saltConcentration);
+    double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0 /*=Rsw*/, saltConcentration);
     rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
     return rho;
 }
@@ -1583,8 +1583,9 @@ applyNumericalAquifers_(const GridView& gridView,
         const auto search = num_aqu_cells.find(cartIx);
         if (search != num_aqu_cells.end()) {
             // numerical aquifer cells are filled with water initially
-            // for co2store the oilphase is used for brine
+            // for co2store the oilphase may be used for the brine
-            const auto watPos = co2store? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx;
+            bool co2store_oil_as_brine = co2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
+            const auto watPos =  co2store_oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx;
             if (FluidSystem::phaseIsActive(watPos)) {
                 this->sat_[watPos][elemIdx] = 1.;
             } else {

flow/flow_ebos_gaswater_dissolution.cpp

@@ -0,0 +1,102 @@
+/*
+  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 3 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/>.
+*/
+#include "config.h"
+
+// Define making clear that the simulator supports AMG
+#define FLOW_SUPPORT_AMG 1
+
+#include <flow/flow_ebos_gaswater_dissolution.hpp>
+
+#include <opm/material/common/ResetLocale.hpp>
+#include <opm/models/blackoil/blackoiltwophaseindices.hh>
+
+#include <opm/grid/CpGrid.hpp>
+#include <opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp>
+#include <opm/simulators/flow/Main.hpp>
+
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
+#include <opm/models/discretization/common/tpfalinearizer.hh>
+
+namespace Opm {
+namespace Properties {
+namespace TTag {
+struct EclFlowGasWaterDissolutionProblem {
+    using InheritsFrom = std::tuple<EclFlowProblem>;
+};
+}
+
+template<class TypeTag>
+struct Linearizer<TypeTag, TTag::EclFlowGasWaterDissolutionProblem> { using type = TpfaLinearizer<TypeTag>; };
+
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::EclFlowGasWaterDissolutionProblem> { using type = BlackOilLocalResidualTPFA<TypeTag>; };
+
+template<class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::EclFlowGasWaterDissolutionProblem> { static constexpr bool value = false; };
+
+template<class TypeTag>
+struct EnableDisgasInWater<TypeTag, TTag::EclFlowGasWaterDissolutionProblem> {
+    static constexpr bool value = true;
+};
+
+//! The indices required by the model
+template<class TypeTag>
+struct Indices<TypeTag, TTag::EclFlowGasWaterDissolutionProblem>
+{
+private:
+    // it is unfortunately not possible to simply use 'TypeTag' here because this leads
+    // to cyclic definitions of some properties. if this happens the compiler error
+    // messages unfortunately are *really* confusing and not really helpful.
+    using BaseTypeTag = TTag::EclFlowProblem;
+    using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
+
+public:
+    typedef BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                    getPropValue<TypeTag, Properties::EnableExtbo>(),
+                                    getPropValue<TypeTag, Properties::EnablePolymer>(),
+                                    getPropValue<TypeTag, Properties::EnableEnergy>(),
+                                    getPropValue<TypeTag, Properties::EnableFoam>(),
+                                    getPropValue<TypeTag, Properties::EnableBrine>(),
+                                    /*PVOffset=*/0,
+                                    /*disabledCompIdx=*/FluidSystem::oilCompIdx,
+                                    getPropValue<TypeTag, Properties::EnableMICP>()> type;
+};
+}}
+
+namespace Opm {
+
+
+// ----------------- Main program -----------------
+int flowEbosGasWaterDissolutionMain(int argc, char** argv, bool outputCout, bool outputFiles)
+{
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    resetLocale();
+
+    FlowMainEbos<Properties::TTag::EclFlowGasWaterDissolutionProblem>
+        mainfunc {argc, argv, outputCout, outputFiles} ;
+    return mainfunc.execute();
+}
+
+int flowEbosGasWaterDissolutionMainStandalone(int argc, char** argv)
+{
+    using TypeTag = Properties::TTag::EclFlowGasWaterDissolutionProblem;
+    auto mainObject = Opm::Main(argc, argv);
+    return mainObject.runStatic<TypeTag>();
+}
+
+}

flow/flow_ebos_gaswater_dissolution.hpp

@@ -0,0 +1,30 @@
+/*
+  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 3 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/>.
+*/
+#ifndef FLOW_EBOS_GASWATER_DISSOLUTION_HPP
+#define FLOW_EBOS_GASWATER_DISSOLUTION_HPP
+
+namespace Opm {
+
+//! \brief Main function used in flow binary.
+int flowEbosGasWaterDissolutionMain(int argc, char** argv, bool outputCout, bool outputFiles);
+
+//! \brief Main function used in flow_gaswater_dissolution binary.
+int flowEbosGasWaterDissolutionMainStandalone(int argc, char** argv);
+
+}
+
+#endif // FLOW_EBOS_GASWATER_DISSOLUTION_HPP

flow/flow_ebos_gaswater_dissolution_diffuse.cpp

@@ -0,0 +1,102 @@
+/*
+  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 3 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/>.
+*/
+#include "config.h"
+
+// Define making clear that the simulator supports AMG
+#define FLOW_SUPPORT_AMG 1
+
+#include <flow/flow_ebos_gaswater_dissolution_diffuse.hpp>
+
+#include <opm/material/common/ResetLocale.hpp>
+#include <opm/models/blackoil/blackoiltwophaseindices.hh>
+
+#include <opm/grid/CpGrid.hpp>
+#include <opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp>
+#include <opm/simulators/flow/Main.hpp>
+
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
+#include <opm/models/discretization/common/tpfalinearizer.hh>
+
+namespace Opm {
+namespace Properties {
+namespace TTag {
+struct EclFlowGasWaterDissolutionDiffuseProblem {
+    using InheritsFrom = std::tuple<EclFlowProblem>;
+};
+}
+
+//template<class TypeTag>
+//struct Linearizer<TypeTag, TTag::EclFlowGasWaterDissolutionProblem> { using type = TpfaLinearizer<TypeTag>; };
+
+//template<class TypeTag>
+//struct LocalResidual<TypeTag, TTag::EclFlowGasWaterDissolutionProblem> { using type = BlackOilLocalResidualTPFA<TypeTag>; };
+
+template<class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::EclFlowGasWaterDissolutionDiffuseProblem> { static constexpr bool value = true; };
+
+template<class TypeTag>
+struct EnableDisgasInWater<TypeTag, TTag::EclFlowGasWaterDissolutionDiffuseProblem> {
+    static constexpr bool value = true;
+};
+
+//! The indices required by the model
+template<class TypeTag>
+struct Indices<TypeTag, TTag::EclFlowGasWaterDissolutionDiffuseProblem>
+{
+private:
+    // it is unfortunately not possible to simply use 'TypeTag' here because this leads
+    // to cyclic definitions of some properties. if this happens the compiler error
+    // messages unfortunately are *really* confusing and not really helpful.
+    using BaseTypeTag = TTag::EclFlowProblem;
+    using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
+
+public:
+    typedef BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                    getPropValue<TypeTag, Properties::EnableExtbo>(),
+                                    getPropValue<TypeTag, Properties::EnablePolymer>(),
+                                    getPropValue<TypeTag, Properties::EnableEnergy>(),
+                                    getPropValue<TypeTag, Properties::EnableFoam>(),
+                                    getPropValue<TypeTag, Properties::EnableBrine>(),
+                                    /*PVOffset=*/0,
+                                    /*disabledCompIdx=*/FluidSystem::oilCompIdx,
+                                    getPropValue<TypeTag, Properties::EnableMICP>()> type;
+};
+}}
+
+namespace Opm {
+
+
+// ----------------- Main program -----------------
+int flowEbosGasWaterDissolutionDiffuseMain(int argc, char** argv, bool outputCout, bool outputFiles)
+{
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    resetLocale();
+
+    FlowMainEbos<Properties::TTag::EclFlowGasWaterDissolutionDiffuseProblem>
+        mainfunc {argc, argv, outputCout, outputFiles} ;
+    return mainfunc.execute();
+}
+
+int flowEbosGasWaterDissolutionDiffuseMainStandalone(int argc, char** argv)
+{
+    using TypeTag = Properties::TTag::EclFlowGasWaterDissolutionDiffuseProblem;
+    auto mainObject = Opm::Main(argc, argv);
+    return mainObject.runStatic<TypeTag>();
+}
+
+}

flow/flow_ebos_gaswater_dissolution_diffuse.hpp

@@ -0,0 +1,30 @@
+/*
+  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 3 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/>.
+*/
+#ifndef FLOW_EBOS_GASWATER_DISSOLUTION_DIFFUSE_HPP
+#define FLOW_EBOS_GASWATER_DISSOLUTION_DIFFUSE_HPP
+
+namespace Opm {
+
+//! \brief Main function used in flow binary.
+int flowEbosGasWaterDissolutionDiffuseMain(int argc, char** argv, bool outputCout, bool outputFiles);
+
+//! \brief Main function used in flow_gaswater_dissolution binary.
+int flowEbosGasWaterDissolutionDiffuseMainStandalone(int argc, char** argv);
+
+}
+
+#endif // FLOW_EBOS_GASWATER_DISSOLUTION_DIFFUSE_HPP

flow/flow_gaswater_dissolution.cpp

@@ -0,0 +1,24 @@
+/*
+  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 3 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/>.
+*/
+#include "config.h"
+#include <flow/flow_ebos_gaswater_dissolution.hpp>
+
+
+int main(int argc, char** argv)
+{
+    return Opm::flowEbosGasWaterDissolutionMainStandalone(argc, argv);
+}

flow/flow_gaswater_dissolution_diffuse.cpp

@@ -0,0 +1,24 @@
+/*
+  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 3 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/>.
+*/
+#include "config.h"
+#include <flow/flow_ebos_gaswater_dissolution_diffuse.hpp>
+
+
+int main(int argc, char** argv)
+{
+    return Opm::flowEbosGasWaterDissolutionDiffuseMainStandalone(argc, argv);
+}

opm/simulators/aquifers/AquiferAnalytical.hpp

@@ -62,6 +62,8 @@ public:
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
     enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
+    enum { has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
+
     enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
 
     static constexpr int numEq = BlackoilIndices::numEq;
@@ -77,6 +79,7 @@ public:
                                           enableEvaporation,
                                           enableBrine,
                                           enableSaltPrecipitation,
+                                          has_disgas_in_water,
                                           BlackoilIndices::numPhases>;
 
     // Constructor

opm/simulators/aquifers/AquiferInterface.hpp

@@ -79,7 +79,8 @@ protected:
 
     int phaseIdx_() const
     {
-        if (co2store_())
+        // If OIL is used to model brine the aquifer should do the same
+        if (co2store_() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
             return FluidSystem::oilPhaseIdx;
 
         return FluidSystem::waterPhaseIdx;

opm/simulators/flow/Main.hpp

@@ -42,6 +42,8 @@
 #include <flow/flow_ebos_onephase_energy.hpp>
 #include <flow/flow_ebos_oilwater_brine.hpp>
 #include <flow/flow_ebos_gaswater_brine.hpp>
+#include <flow/flow_ebos_gaswater_dissolution.hpp>
+#include <flow/flow_ebos_gaswater_dissolution_diffuse.hpp>
 #include <flow/flow_ebos_energy.hpp>
 #include <flow/flow_ebos_oilwater_polymer.hpp>
 #include <flow/flow_ebos_oilwater_polymer_injectivity.hpp>
@@ -632,6 +634,7 @@ private:
     int runTwoPhase(const Phases& phases)
     {
         const bool diffusive = eclipseState_->getSimulationConfig().isDiffusive();
+        const bool disgasw = eclipseState_->getSimulationConfig().hasDISGASW();
 
         // oil-gas
         if (phases.active( Phase::OIL ) && phases.active( Phase::GAS )) {
@@ -655,12 +658,19 @@ private:
 
         // gas-water
         else if ( phases.active( Phase::GAS ) && phases.active( Phase::WATER ) ) {
+            if (disgasw) {
+                if (diffusive) {
+                    return flowEbosGasWaterDissolutionDiffuseMain(argc_, argv_, outputCout_, outputFiles_);
+                }
+                return flowEbosGasWaterDissolutionMain(argc_, argv_, outputCout_, outputFiles_);
+            }
             if (diffusive) {
                 if (outputCout_) {
-                    std::cerr << "The DIFFUSE option is not available for the two-phase gas/water model." << std::endl;
+                    std::cerr << "The DIFFUSE option is not available for the two-phase gas/water model without disgasw." << std::endl;
                 }
                 return EXIT_FAILURE;
             }
+
             return flowEbosGasWaterMain(argc_, argv_, outputCout_, outputFiles_);
         }
         else {

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -305,6 +305,8 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
                 (pvtReg, segment_pressure[seg],
                  std::max(0.0, Rs),
                  std::max(0.0, Rv),
+                 0.0, // Rsw
+                 0.0, // Rvw
                  temperature, saltConc, surf_rates, resv_rates);
         }
 
@@ -340,7 +342,7 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
 
         if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
             segments.phase_viscosity[seg * well_.numPhases() + pu.phase_pos[Water]] =
-                FluidSystem::waterPvt().viscosity(pvtReg, temperature, segment_pressure[seg], saltConc);
+                FluidSystem::waterPvt().viscosity(pvtReg, temperature, segment_pressure[seg],  0.0 /*Rsw*/, saltConc);
         }
 
         if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {

opm/simulators/wells/MultisegmentWellSegments.cpp

@@ -143,11 +143,12 @@ computeFluidProperties(const EvalWell& temperature,
 
         const EvalWell seg_pressure = primary_variables.getSegmentPressure(seg);
         if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            EvalWell rsw(0.0);
             const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
             b[waterCompIdx] =
-                FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration);
+                FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rsw, saltConcentration);
             visc[waterCompIdx] =
-                FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure, saltConcentration);
+                FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure, rsw, saltConcentration);
             // TODO: double check here
             // TODO: should not we use phaseIndex here?
             phase_densities[waterCompIdx] = b[waterCompIdx] * surf_dens[waterCompIdx];
@@ -348,10 +349,12 @@ getSurfaceVolume(const EvalWell& temperature,
     std::vector<EvalWell> b(well_.numComponents(), 0.);
     if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
         const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
+        EvalWell rsw(0.0);
         b[waterCompIdx] =
             FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index,
                                                                  temperature,
                                                                  seg_pressure,
+                                                                 rsw,
                                                                  saltConcentration);
     }
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -84,6 +84,10 @@ namespace Opm
             OPM_THROW(std::runtime_error, "dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet."
             << " \n See  (WCONINJE item 10 / WCONHIST item 8)");
         }
+        if constexpr (!Indices::oilEnabled && Indices::numPhases > 1) {
+            OPM_THROW(std::runtime_error, "water + gas case not supported by multisegment well yet");
+        }
+
     }
 
 
@@ -1546,6 +1550,7 @@ namespace Opm
 
         auto& ws = well_state.well(this->index_of_well_);
         ws.dissolved_gas_rate = 0;
+        ws.dissolved_gas_rate_in_water = 0;
         ws.vaporized_oil_rate = 0;
         ws.vaporized_wat_rate = 0;
 

opm/simulators/wells/RateConverter.hpp

@@ -104,6 +104,8 @@ namespace Opm {
                     ra.rv = 0.0;
                     ra.pv = 0.0;
                     ra.saltConcentration = 0.0;
+                    ra.rsw = 0.0;
+                    ra.rvw = 0.0;
 
                 }
 
@@ -162,6 +164,12 @@ namespace Opm {
                             attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
                         }
                         attr.saltConcentration += fs.saltConcentration().value() * hydrocarbonPV;
+                        if (FluidSystem::enableDissolvedGasInWater()) {
+                            attr.rsw += fs.Rsw().value() * hydrocarbonPV; // scale with total volume?
+                        }
+                        if (FluidSystem::enableVaporizedWater()) {
+                            attr.rvw += fs.Rvw().value() * hydrocarbonPV; // scale with total volume?
+                        }
                     }
 
                     if (pv_cell > 0.) {
@@ -183,6 +191,12 @@ namespace Opm {
                             attr.temperature += fs.temperature(FluidSystem::waterPhaseIdx).value() * pv_cell;
                         }
                         attr.saltConcentration += fs.saltConcentration().value() * pv_cell;
+                        if (FluidSystem::enableDissolvedGasInWater()) {
+                            attr.rsw += fs.Rsw().value() * pv_cell; 
+                        }
+                        if (FluidSystem::enableVaporizedWater()) {
+                            attr.rvw += fs.Rvw().value() * pv_cell;
+                        }
                     }
                 }
 
@@ -199,11 +213,15 @@ namespace Opm {
                           const double rs_hpv_sum = comm.sum(attri_hpv.rs);
                           const double rv_hpv_sum = comm.sum(attri_hpv.rv);
                           const double sc_hpv_sum = comm.sum(attri_hpv.saltConcentration);
+                          const double rsw_hpv_sum = comm.sum(attri_hpv.rsw);
+                          const double rvw_hpv_sum = comm.sum(attri_hpv.rvw);
 
                           ra.pressure = p_hpv_sum / hpv_sum;
                           ra.temperature = T_hpv_sum / hpv_sum;
                           ra.rs = rs_hpv_sum / hpv_sum;
                           ra.rv = rv_hpv_sum / hpv_sum;
+                          ra.rsw = rsw_hpv_sum / hpv_sum;
+                          ra.rvw = rvw_hpv_sum / hpv_sum;
                           ra.pv = hpv_sum;
                           ra.saltConcentration = sc_hpv_sum / hpv_sum;
                       } else {
@@ -215,12 +233,16 @@ namespace Opm {
                           const double T_pv_sum = comm.sum(attri_pv.temperature);
                           const double rs_pv_sum = comm.sum(attri_pv.rs);
                           const double rv_pv_sum = comm.sum(attri_pv.rv);
+                          const double rsw_pv_sum = comm.sum(attri_pv.rsw);
+                          const double rvw_pv_sum = comm.sum(attri_pv.rvw);
                           const double sc_pv_sum = comm.sum(attri_pv.saltConcentration);
 
                           ra.pressure = p_pv_sum / pv_sum;
                           ra.temperature = T_pv_sum / pv_sum;
                           ra.rs = rs_pv_sum / pv_sum;
                           ra.rv = rv_pv_sum / pv_sum;
+                          ra.rsw = rsw_pv_sum / pv_sum;
+                          ra.rvw = rvw_pv_sum / pv_sum;
                           ra.pv = pv_sum;
                           ra.saltConcentration = sc_pv_sum / pv_sum;
                       }
@@ -277,12 +299,24 @@ namespace Opm {
 
                 std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
 
+                // Actual Rsw and Rvw:
+                double Rsw = ra.rsw;
+                double Rvw = ra.rvw;
+                // Determinant of 'R' matrix
+                const double detRw = 1.0 - (Rsw * Rvw);
+
                 if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
-                    // q[w]_r = q[w]_s / bw
+                    // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s)
 
-                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
+                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rsw, saltConcentration);
 
-                    coeff[iw] = 1.0 / bw;
+                    const double den = bw * detRw;
+
+                    coeff[iw] += 1.0 / den;
+
+                    if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
+                        coeff[ig] -= ra.rvw / den;
+                    }
                 }
 
                 // Actual Rs and Rv:
@@ -292,6 +326,13 @@ namespace Opm {
                 // Determinant of 'R' matrix
                 const double detR = 1.0 - (Rs * Rv);
 
+                // Currently we only support either gas in water or gas in oil
+                // not both
+                if(detR != 1 && detRw != 1 ) {
+                    std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR);
+                    throw(msg);
+                }
+
                 if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
                     // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
 
@@ -307,14 +348,19 @@ namespace Opm {
 
                 if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
                     // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
-
+                    const double bg  = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv, Rvw);
-                    const double bg  = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv, 0.0 /*=Rvw*/);
+                    if (FluidSystem::enableDissolvedGasInWater()) {
-                    const double den = bg * detR;
+                        const double denw = bg * detRw;
-
+                        coeff[ig] += 1.0 / denw;
-                    coeff[ig] += 1.0 / den;
+                        if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
-
+                           coeff[iw] -= ra.rsw / denw;
-                    if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
+                        }
-                        coeff[io] -= ra.rs / den;
+                    } else {
+                        const double den = bg * detR;
+                        coeff[ig] += 1.0 / den;
+                        if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
+                            coeff[io] -= ra.rs / den;
+                        }
                     }
                 }
             }
@@ -338,7 +384,7 @@ namespace Opm {
                 if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
                     // q[w]_r = q[w]_s / bw
 
-                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
+                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, saltConcentration);
 
                     coeff[iw] = 1.0 / bw;
                 }
@@ -385,7 +431,8 @@ namespace Opm {
                 const auto& ra = this->attr_.attributes(r);
 
                 this->calcReservoirVoidageRates(pvtRegionIdx,
-                                                ra.pressure, ra.rs, ra.rv,
+                                                ra.pressure, ra.rs, ra.rv, 
+                                                ra.rsw, ra.rvw, 
                                                 ra.temperature,
                                                 ra.saltConcentration,
                                                 surface_rates,
@@ -409,6 +456,10 @@ namespace Opm {
              * \param[in] rs Dissolved gas/oil ratio.
              *
              * \param[in] rv Vaporised oil/gas ratio.
+             *        
+             * \param[in] rsw Dissolved gas/water ratio.
+             *
+             * \param[in] rwv Vaporised water/gas ratio.
              *
              * \param[in] T Temperature.  Unused in non-thermal simulation
              *    runs.
@@ -427,6 +478,8 @@ namespace Opm {
                                            const double        p,
                                            const double        rs,
                                            const double        rv,
+                                           const double        rsw,
+                                           const double        rvw,
                                            const double        T,
                                            const double        saltConcentration,
                                            const SurfaceRates& surface_rates,
@@ -440,15 +493,29 @@ namespace Opm {
                 const auto [Rs, Rv] = this->
                     dissolvedVaporisedRatio(io, ig, rs, rv, surface_rates);
 
+                const auto [Rsw, Rvw] = this->
+                    dissolvedVaporisedRatio(iw, ig, rsw, rvw, surface_rates);
+
+
                 std::fill_n(&voidage_rates[0], pu.num_phases, 0.0);
 
+
+                // Determinant of 'R' matrix
+                const auto detRw = 1.0 - (Rsw * Rvw);
+
                 if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
-                    // q[w]_r = q[w]_s / bw
+                    // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s)
+                    voidage_rates[iw] = surface_rates[iw];
+
                     const auto bw = FluidSystem::waterPvt()
                         .inverseFormationVolumeFactor(pvtRegionIdx, T, p,
+                                                      Rsw,
                                                       saltConcentration);
 
-                    voidage_rates[iw] = surface_rates[iw] / bw;
+                    if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
+                        voidage_rates[iw] -= Rvw * surface_rates[ig];
+                    }
+                    voidage_rates[iw] /= bw * detRw;
                 }
 
                 // Determinant of 'R' matrix
@@ -467,18 +534,32 @@ namespace Opm {
                     voidage_rates[io] /= bo * detR;
                 }
 
+                // we only support either gas in water
+                // or gas in oil
+                if(detR != 1 && detRw != 1 ) {
+                    std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR);
+                    throw(msg);
+                }
                 if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
                     // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
                     voidage_rates[ig] = surface_rates[ig];
                     if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
                         voidage_rates[ig] -= Rs * surface_rates[io];
                     }
+                    if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
+                        voidage_rates[ig] -= Rsw * surface_rates[iw];
+                    }
 
                     const auto bg = FluidSystem::gasPvt()
                         .inverseFormationVolumeFactor(pvtRegionIdx, T, p,
-                                                      Rv, 0.0 /*=Rvw*/);
+                                                      Rv, Rvw);
 
-                    voidage_rates[ig] /= bg * detR;
+                    // we only support either gas in water or gas in oil
+                    if (detRw == 1) {
+                        voidage_rates[ig] /= bg * detR;
+                    } else { 
+                        voidage_rates[ig] /= bg * detRw;
+                    }
                 }
             }
 
@@ -537,6 +618,8 @@ namespace Opm {
                     , temperature(0.0)
                     , rs(0.0)
                     , rv(0.0)
+                    , rsw(0.0)
+                    , rvw(0.0)
                     , pv(0.0)
                     , saltConcentration(0.0)
                 {}
@@ -545,6 +628,8 @@ namespace Opm {
                 double temperature;
                 double rs;
                 double rv;
+                double rsw;
+                double rvw;
                 double pv;
                 double saltConcentration;
             };

opm/simulators/wells/SingleWellState.hpp

@@ -55,6 +55,7 @@ public:
     double thp{0};
     double temperature{0};
     double dissolved_gas_rate{0};
+    double dissolved_gas_rate_in_water{0};
     double vaporized_oil_rate{0};
     double vaporized_wat_rate{0};
     std::vector<double> well_potentials;

opm/simulators/wells/StandardWell.hpp

@@ -266,6 +266,7 @@ namespace Opm
                                                          std::vector<double>& rsmax_perf,
                                                          std::vector<double>& rvmax_perf,
                                                          std::vector<double>& rvwmax_perf,
+                                                         std::vector<double>& rswmax_perf,
                                                          std::vector<double>& surf_dens_perf) const;
 
         void computeWellConnectionDensitesPressures(const Simulator& ebosSimulator,
@@ -274,6 +275,7 @@ namespace Opm
                                                     const std::vector<double>& rsmax_perf,
                                                     const std::vector<double>& rvmax_perf,
                                                     const std::vector<double>& rvwmax_perf,
+                                                    const std::vector<double>& rswmax_perf,
                                                     const std::vector<double>& surf_dens_perf,
                                                     DeferredLogger& deferred_logger);
 
@@ -289,6 +291,7 @@ namespace Opm
                                  const bool allow_cf,
                                  std::vector<EvalWell>& cq_s,
                                  double& perf_dis_gas_rate,
+                                 double& perf_dis_gas_rate_in_water,
                                  double& perf_vap_oil_rate,
                                  double& perf_vap_wat_rate,
                                  DeferredLogger& deferred_logger) const;
@@ -309,6 +312,7 @@ namespace Opm
                              const Value& rs,
                              const Value& rv,
                              const Value& rvw,
+                             const Value& rsw,
                              std::vector<Value>& b_perfcells_dense,
                              const double Tw,
                              const int perf,
@@ -317,6 +321,7 @@ namespace Opm
                              const std::vector<Value>& cmix_s,
                              std::vector<Value>& cq_s,
                              double& perf_dis_gas_rate,
+                             double& perf_dis_gas_rate_in_water,
                              double& perf_vap_oil_rate,
                              double& perf_vap_wat_rate,
                              DeferredLogger& deferred_logger) const;

opm/simulators/wells/StandardWellConnections.cpp

@@ -92,6 +92,7 @@ computeDensities(const std::vector<Scalar>& perfComponentRates,
                  const std::vector<Scalar>& rsmax_perf,
                  const std::vector<Scalar>& rvmax_perf,
                  const std::vector<Scalar>& rvwmax_perf,
+                 const std::vector<Scalar>& rswmax_perf,
                  const std::vector<Scalar>& surf_dens_perf,
                  DeferredLogger& deferred_logger)
 {
@@ -216,7 +217,7 @@ computeDensities(const std::vector<Scalar>& perfComponentRates,
             if (d <= 0.0) {
                 std::ostringstream sstr;
                 sstr << "Problematic d value " << d << " obtained for well " << well_.name()
-                     << " during ccomputeConnectionDensities with rs " << rs
+                     << " during computeConnectionDensities with rs " << rs
                      << ", rv " << rv
                      << " obtaining d " << d
                      << " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
@@ -232,34 +233,44 @@ computeDensities(const std::vector<Scalar>& perfComponentRates,
                     x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/d;
                 }
             }
-            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+        }
-                //matrix system: (mix[oilpos] = q_os, x[oilpos] = bo*q_or, etc...)
+
-                //┌             ┐   ┌                ┐  ┌           ┐
+        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                //│mix[oilpos]  │   | 1     Rv     0 |  |x[oilpos]  |
+            //matrix system: (mix[oilpos] = q_os, x[oilpos] = bo*q_or, etc...)
-                //│mix[gaspos]  │ = │ Rs    1      0 │  │x[gaspos]  │
+            //┌             ┐   ┌                ┐  ┌           ┐
-                //│mix[waterpos]│   │ 0     Rvw    1 │  │x[waterpos │
+            //│mix[oilpos]  │   | 1     Rv     0 |  |x[oilpos]  |
-                //└             ┘   └                ┘  └           ┘
+            //│mix[gaspos]  │ = │ Rs    1     Rsw│  │x[gaspos]  │
-                const unsigned waterpos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
+            //│mix[waterpos]│   │ 0     Rvw    1 │  │x[waterpos │
-                Scalar rvw = 0.0;
+            //└             ┘   └                ┘  └           ┘
-                if (!rvwmax_perf.empty() && mix[gaspos] > 1e-12) {
-                    rvw = std::min(mix[waterpos]/mix[gaspos], rvwmax_perf[perf]);
-                }
-                if (rvw > 0.0) {
-                    // Subtract water in gas from water mixture
-                    x[waterpos] = mix[waterpos] - x[gaspos] * rvw;
-                }
-            }
-        } else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-            //no oil
-            const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
             const unsigned waterpos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
+            const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
             Scalar rvw = 0.0;
             if (!rvwmax_perf.empty() && mix[gaspos] > 1e-12) {
                 rvw = std::min(mix[waterpos]/mix[gaspos], rvwmax_perf[perf]);
             }
-            if (rvw > 0.0) {
+            Scalar rsw = 0.0;
-               // Subtract water in gas from water mixture
+            if (!rswmax_perf.empty() && mix[waterpos] > 1e-12) {
-               x[waterpos] = mix[waterpos] - mix[gaspos] * rvw;
+                rsw = std::min(mix[gaspos]/mix[waterpos], rswmax_perf[perf]);
+            }
+            const Scalar d = 1.0 - rsw*rvw;
+            if (d <= 0.0) {
+                std::ostringstream sstr;
+                sstr << "Problematic d value " << d << " obtained for well " << well_.name()
+                    << " during computeConnectionDensities with rsw " << rsw
+                    << ", rvw " << rvw
+                    << " obtaining d " << d
+                    << " Continue as if no dissolution (rsw = 0) and vaporization (rvw = 0) "
+                    << " for this connection.";
+                deferred_logger.debug(sstr.str());
+            } else {
+                if (rsw > 0.0) {
+                    // Subtract gas in water from gas mixture
+                    x[gaspos] = (mix[gaspos] - mix[waterpos]*rsw)/d;
+                }
+                if (rvw > 0.0) {
+                    // Subtract water in gas from water mixture
+                    x[waterpos] = (mix[waterpos] - mix[gaspos]*rvw)/d;
+                }
             }
         }
 
@@ -288,6 +299,7 @@ computePropertiesForPressures(const WellState& well_state,
                               std::vector<Scalar>& rsmax_perf,
                               std::vector<Scalar>& rvmax_perf,
                               std::vector<Scalar>& rvwmax_perf,
+                              std::vector<Scalar>& rswmax_perf,
                               std::vector<Scalar>& surf_dens_perf) const
 {
     const int nperf = well_.numPerfs();
@@ -312,6 +324,7 @@ computePropertiesForPressures(const WellState& well_state,
       //rvw is only used if both water and gas is present
     if (waterPresent && gasPresent) {
         rvwmax_perf.resize(nperf);
+        rswmax_perf.resize(nperf);
     }
 
     // Compute the average pressure in each well block
@@ -330,8 +343,21 @@ computePropertiesForPressures(const WellState& well_state,
 
         if (waterPresent) {
             const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
-            b_perf[ waterCompIdx + perf * well_.numComponents()] =
+            double rsw = 0.0;
-            FluidSystem::waterPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, saltConcentration);
+            if (FluidSystem::enableDissolvedGasInWater()) {
+                // TODO support mutual solubility in water and oil
+                assert(!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
+                const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]);
+                rswmax_perf[perf] = FluidSystem::waterPvt().saturatedGasDissolutionFactor(region_idx, temperature, p_avg, saltConcentration);
+                if (waterrate > 0) {
+                    const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0);
+                    if (gasrate > 0) {
+                        rsw = waterrate / gasrate;
+                    }
+                    rsw = std::min(rsw, rswmax_perf[perf]);
+                }
+            }
+            b_perf[ waterCompIdx + perf * well_.numComponents()] = FluidSystem::waterPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, rsw, saltConcentration);
         }
 
         if (gasPresent) {
@@ -456,6 +482,7 @@ computeProperties(const WellState& well_state,
                   const std::vector<Scalar>& rsmax_perf,
                   const std::vector<Scalar>& rvmax_perf,
                   const std::vector<Scalar>& rvwmax_perf,
+                  const std::vector<Scalar>& rswmax_perf,
                   const std::vector<Scalar>& surf_dens_perf,
                   DeferredLogger& deferred_logger)
 {
@@ -521,7 +548,7 @@ computeProperties(const WellState& well_state,
         }
     }
 
-    this->computeDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger);
+    this->computeDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, rswmax_perf, surf_dens_perf, deferred_logger);
     this->computePressureDelta();
 }
 

opm/simulators/wells/StandardWellConnections.hpp

@@ -49,6 +49,7 @@ public:
                                        std::vector<Scalar>& rsmax_perf,
                                        std::vector<Scalar>& rvmax_perf,
                                        std::vector<Scalar>& rvwmax_perf,
+                                       std::vector<Scalar>& rswmax_perf,
                                        std::vector<Scalar>& surf_dens_perf) const;
 
     //! \brief Compute connection properties (densities, pressure drop, ...)
@@ -61,6 +62,7 @@ public:
                            const std::vector<Scalar>& rsmax_perf,
                            const std::vector<Scalar>& rvmax_perf,
                            const std::vector<Scalar>& rvwmax_perf,
+                           const std::vector<Scalar>& rswmax_perf,
                            const std::vector<Scalar>& surf_dens_perf,
                            DeferredLogger& deferred_logger);
 
@@ -84,6 +86,7 @@ private:
                           const std::vector<Scalar>& rsmax_perf,
                           const std::vector<Scalar>& rvmax_perf,
                           const std::vector<Scalar>& rvwmax_perf,
+                          const std::vector<Scalar>& rswmax_perf,
                           const std::vector<Scalar>& surf_dens_perf,
                           DeferredLogger& deferred_logger);
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -97,6 +97,7 @@ namespace Opm
                         const bool allow_cf,
                         std::vector<EvalWell>& cq_s,
                         double& perf_dis_gas_rate,
+                        double& perf_dis_gas_rate_in_water,
                         double& perf_vap_oil_rate,
                         double& perf_vap_wat_rate,
                         DeferredLogger& deferred_logger) const
@@ -106,6 +107,8 @@ namespace Opm
         const EvalWell rs = this->extendEval(fs.Rs());
         const EvalWell rv = this->extendEval(fs.Rv());
         const EvalWell rvw = this->extendEval(fs.Rvw());
+        const EvalWell rsw = this->extendEval(fs.Rsw());
+
 
         std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
         for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
@@ -147,6 +150,7 @@ namespace Opm
                         rs,
                         rv,
                         rvw,
+                        rsw,
                         b_perfcells_dense,
                         Tw,
                         perf,
@@ -155,6 +159,7 @@ namespace Opm
                         cmix_s,
                         cq_s,
                         perf_dis_gas_rate,
+                        perf_dis_gas_rate_in_water,
                         perf_vap_oil_rate,
                         perf_vap_wat_rate,
                         deferred_logger);
@@ -177,6 +182,7 @@ namespace Opm
         const Scalar rs = fs.Rs().value();
         const Scalar rv = fs.Rv().value();
         const Scalar rvw = fs.Rvw().value();
+        const Scalar rsw = fs.Rsw().value();
         std::vector<Scalar> b_perfcells_dense(this->num_components_, 0.0);
         for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx)) {
@@ -208,6 +214,7 @@ namespace Opm
         Scalar perf_dis_gas_rate = 0.0;
         Scalar perf_vap_oil_rate = 0.0;
         Scalar perf_vap_wat_rate = 0.0;
+        Scalar perf_dis_gas_rate_in_water = 0.0;
 
         // surface volume fraction of fluids within wellbore
         std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
@@ -221,6 +228,7 @@ namespace Opm
                         rs,
                         rv,
                         rvw,
+                        rsw,
                         b_perfcells_dense,
                         Tw,
                         perf,
@@ -229,6 +237,7 @@ namespace Opm
                         cmix_s,
                         cq_s,
                         perf_dis_gas_rate,
+                        perf_dis_gas_rate_in_water,
                         perf_vap_oil_rate,
                         perf_vap_wat_rate,
                         deferred_logger);
@@ -244,6 +253,7 @@ namespace Opm
                     const Value& rs,
                     const Value& rv,
                     const Value& rvw,
+                    const Value& rsw,
                     std::vector<Value>& b_perfcells_dense,
                     const double Tw,
                     const int perf,
@@ -252,6 +262,7 @@ namespace Opm
                     const std::vector<Value>& cmix_s,
                     std::vector<Value>& cq_s,
                     double& perf_dis_gas_rate,
+                    double& perf_dis_gas_rate_in_water,
                     double& perf_vap_oil_rate,
                     double& perf_vap_wat_rate,
                     DeferredLogger& deferred_logger) const
@@ -292,12 +303,20 @@ namespace Opm
                     perf_dis_gas_rate = getValue(dis_gas);
                     perf_vap_oil_rate = getValue(vap_oil);
                 }
-                if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            }
-                    const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
+
-                    const Value vap_wat = rvw * cq_sGas;
+            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                    cq_s[waterCompIdx] += vap_wat;
+                const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
-                    if (this->isProducer())
+                const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
-                        perf_vap_wat_rate = getValue(vap_wat);
+                const Value cq_sWat = cq_s[waterCompIdx];
+                const Value cq_sGas = cq_s[gasCompIdx];
+                const Value vap_wat = rvw * cq_sGas;
+                const Value dis_gas_wat = rsw * cq_sWat;
+                cq_s[waterCompIdx] += vap_wat;
+                cq_s[gasCompIdx]   += dis_gas_wat;
+                if (this->isProducer()) {
+                    perf_vap_wat_rate = getValue(vap_wat);
+                    perf_dis_gas_rate_in_water = getValue(dis_gas_wat);
                 }
             }
 
@@ -319,9 +338,32 @@ namespace Opm
             // compute volume ratio between connection at standard conditions
             Value volumeRatio = bhp * 0.0; // initialize it with the correct type
 ;
-            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
                 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
-                volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
+                const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+                // Incorporate RSW/RVW factors if both water and gas active
+                const Value d = 1.0 - rvw * rsw;
+
+                if (d <= 0.0) {
+                    std::ostringstream sstr;
+                    sstr << "Problematic d value " << d << " obtained for well " << this->name()
+                         << " during computePerfRate calculations with rsw " << rsw
+                         << ", rvw " << rvw << " and pressure " << pressure
+                         << " obtaining d " << d
+                         << " Continue as if no dissolution (rsw = 0) and vaporization (rvw = 0) "
+                         << " for this connection.";
+                    deferred_logger.debug(sstr.str());
+                }
+                const Value tmp_wat = d > 0.0? (cmix_s[waterCompIdx] - rvw * cmix_s[gasCompIdx]) / d : cmix_s[waterCompIdx];
+                volumeRatio += tmp_wat / b_perfcells_dense[waterCompIdx];
+
+                const Value tmp_gas =  d > 0.0? (cmix_s[gasCompIdx] - rsw * cmix_s[waterCompIdx]) / d : cmix_s[waterCompIdx];
+                volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
+            } else {
+                if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                    const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
+                    volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
+                }
             }
 
             if constexpr (Indices::enableSolvent) {
@@ -407,10 +449,12 @@ namespace Opm
                         perf_vap_wat_rate = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
                     }
                 }
-                else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                     //no oil
                     const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+                    const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
                     perf_vap_wat_rate = getValue(rvw) * getValue(cq_s[gasCompIdx]);
+                    perf_dis_gas_rate_in_water = getValue(rsw) * getValue(cq_s[waterCompIdx]);
                 }
             }
         }
@@ -458,6 +502,7 @@ namespace Opm
 
         ws.vaporized_oil_rate = 0;
         ws.dissolved_gas_rate = 0;
+        ws.dissolved_gas_rate_in_water = 0;
         ws.vaporized_wat_rate = 0;
 
         const int np = this->number_of_phases_;
@@ -517,6 +562,7 @@ namespace Opm
         {
             const auto& comm = this->parallel_well_info_.communication();
             ws.dissolved_gas_rate = comm.sum(ws.dissolved_gas_rate);
+            ws.dissolved_gas_rate_in_water = comm.sum(ws.dissolved_gas_rate_in_water);
             ws.vaporized_oil_rate = comm.sum(ws.vaporized_oil_rate);
             ws.vaporized_wat_rate = comm.sum(ws.vaporized_wat_rate);
         }
@@ -584,12 +630,13 @@ namespace Opm
         getMobilityEval(ebosSimulator, perf, mob, deferred_logger);
 
         double perf_dis_gas_rate = 0.;
+        double perf_dis_gas_rate_in_water = 0.;
         double perf_vap_oil_rate = 0.;
         double perf_vap_wat_rate = 0.;
         double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants,  cell_idx);
         const double Tw = this->well_index_[perf] * trans_mult;
         computePerfRateEval(intQuants, mob, bhp, Tw, perf, allow_cf,
-                            cq_s, perf_dis_gas_rate, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
+                            cq_s, perf_dis_gas_rate, perf_dis_gas_rate_in_water, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
 
         auto& ws = well_state.well(this->index_of_well_);
         auto& perf_data = ws.perf_data;
@@ -608,6 +655,7 @@ namespace Opm
         // updating the solution gas rate and solution oil rate
         if (this->isProducer()) {
             ws.dissolved_gas_rate += perf_dis_gas_rate;
+            ws.dissolved_gas_rate_in_water += perf_dis_gas_rate_in_water;
             ws.vaporized_oil_rate += perf_vap_oil_rate;
             ws.vaporized_wat_rate += perf_vap_wat_rate;
         }
@@ -1288,6 +1336,7 @@ namespace Opm
                                                 std::vector<double>& rsmax_perf,
                                                 std::vector<double>& rvmax_perf,
                                                 std::vector<double>& rvwmax_perf,
+                                                std::vector<double>& rswmax_perf,
                                                 std::vector<double>& surf_dens_perf) const
     {
         std::function<Scalar(int,int)> getTemperature =
@@ -1326,6 +1375,7 @@ namespace Opm
                                                          rsmax_perf,
                                                          rvmax_perf,
                                                          rvwmax_perf,
+                                                         rswmax_perf,
                                                          surf_dens_perf);
     }
 
@@ -1448,6 +1498,7 @@ namespace Opm
                                            const std::vector<double>& rsmax_perf,
                                            const std::vector<double>& rvmax_perf,
                                            const std::vector<double>& rvwmax_perf,
+                                           const std::vector<double>& rswmax_perf,
                                            const std::vector<double>& surf_dens_perf,
                                            DeferredLogger& deferred_logger)
     {
@@ -1481,6 +1532,7 @@ namespace Opm
                                              rsmax_perf,
                                              rvmax_perf,
                                              rvwmax_perf,
+                                             rswmax_perf,
                                              surf_dens_perf,
                                              deferred_logger);
     }
@@ -1503,9 +1555,10 @@ namespace Opm
          std::vector<double> rsmax_perf;
          std::vector<double> rvmax_perf;
          std::vector<double> rvwmax_perf;
+         std::vector<double> rswmax_perf;
          std::vector<double> surf_dens_perf;
-         computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf);
+         computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, rswmax_perf, surf_dens_perf);
-         computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger);
+         computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, rswmax_perf, surf_dens_perf, deferred_logger);
     }
 
 
@@ -1927,12 +1980,13 @@ namespace Opm
 
             std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
             double perf_dis_gas_rate = 0.;
+            double perf_dis_gas_rate_in_water = 0.;
             double perf_vap_oil_rate = 0.;
             double perf_vap_wat_rate = 0.;
             double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx);
             const double Tw = this->well_index_[perf] * trans_mult;
             computePerfRateEval(int_quant, mob, bhp, Tw, perf, allow_cf,
-                                cq_s, perf_dis_gas_rate, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
+                                cq_s, perf_dis_gas_rate, perf_dis_gas_rate_in_water, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
             // TODO: make area a member
             const double area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
             const auto& material_law_manager = ebos_simulator.problem().materialLawManager();

opm/simulators/wells/WellInterface.hpp

@@ -113,6 +113,7 @@ public:
     static constexpr bool has_foam = getPropValue<TypeTag, Properties::EnableFoam>();
     static constexpr bool has_brine = getPropValue<TypeTag, Properties::EnableBrine>();
     static constexpr bool has_watVapor = getPropValue<TypeTag, Properties::EnableEvaporation>();
+    static constexpr bool has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>();
     static constexpr bool has_saltPrecip = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
     static constexpr bool has_micp = getPropValue<TypeTag, Properties::EnableMICP>();
 
@@ -125,6 +126,7 @@ public:
                                           has_watVapor,
                                           has_brine,
                                           has_saltPrecip,
+                                          has_disgas_in_water,
                                           Indices::numPhases >;
     /// Constructor
     WellInterface(const Well& well,

opm/simulators/wells/WellState.cpp

@@ -495,7 +495,7 @@ WellState::report(const int* globalCellIdxMap,
             well.rates.set(rt::alq, 0.0);
         }
 
-        well.rates.set(rt::dissolved_gas, ws.dissolved_gas_rate);
+        well.rates.set(rt::dissolved_gas, ws.dissolved_gas_rate + ws.dissolved_gas_rate_in_water);
         well.rates.set(rt::vaporized_oil, ws.vaporized_oil_rate);
         well.rates.set(rt::vaporized_water, ws.vaporized_wat_rate);