add output of dissolved solvent rswsol

ebos/eclgenericoutputblackoilmodule.cc

@@ -474,6 +474,7 @@ assignToSolution(data::Solution& sol)
         DataEntry{"PRESPOTF", UnitSystem::measure::pressure,           mechPotentialPressForce_},
         DataEntry{"PRES_OVB", UnitSystem::measure::pressure,           overburdenPressure_},
         DataEntry{"RSW",      UnitSystem::measure::gas_oil_ratio,      rsw_},
+        DataEntry{"RSWSOL",   UnitSystem::measure::gas_oil_ratio,      rswSol_},
         DataEntry{"RVW",      UnitSystem::measure::oil_gas_ratio,      rvw_},
         DataEntry{"SALTP",    UnitSystem::measure::identity,           pSalt_},
         DataEntry{"SS_X",     UnitSystem::measure::identity,           extboX_},
@@ -652,6 +653,12 @@ setRestart(const data::Solution& sol,
         so -= sSol_[elemIdx];
     }
 
+    if (!rswSol_.empty()) {
+        if (sol.has("RSWSOL"))
+            rswSol_[elemIdx] = sol.data<Scalar>("RSWSOL")[globalDofIndex];
+
+    }
+
     assert(!saturation_[oilPhaseIdx].empty());
     saturation_[oilPhaseIdx][elemIdx] = so;
 
@@ -966,6 +973,9 @@ doAllocBuffers(const unsigned bufferSize,
 
     if (enableSolvent_) {
         sSol_.resize(bufferSize, 0.0);
+        if (eclState_.getSimulationConfig().hasDISGASW()) {
+            rswSol_.resize(bufferSize, 0.0);
+        }
     }
 
     if (enablePolymer_) {

ebos/eclgenericoutputblackoilmodule.hh

@@ -106,6 +106,14 @@ public:
         return 0;
     }
 
+    Scalar getSolventRsw(unsigned elemIdx) const
+    {
+        if (rswSol_.size() > elemIdx)
+            return rswSol_[elemIdx];
+
+        return 0;
+    }
+
     Scalar getPolymerConcentration(unsigned elemIdx) const
     {
         if (cPolymer_.size() > elemIdx)
@@ -395,6 +403,7 @@ protected:
     ScalarBuffer oilSaturationPressure_;
     ScalarBuffer drsdtcon_;
     ScalarBuffer sSol_;
+    ScalarBuffer rswSol_;
     ScalarBuffer cPolymer_;
     ScalarBuffer cFoam_;
     ScalarBuffer cSalt_;

ebos/eclgenericproblem.hh

@@ -167,6 +167,11 @@ public:
      */
     Scalar solventSaturation(unsigned elemIdx) const;
 
+    /*!
+     * \brief Returns the initial solvent dissolved in water for a given a cell index
+     */
+    Scalar solventRsw(unsigned elemIdx) const;
+
     /*!
      * \brief Returns the dynamic drsdt convective mixing value
      */
@@ -280,6 +285,7 @@ public:
         serializer(minOilPressure_);
         serializer(overburdenPressure_);
         serializer(solventSaturation_);
+        serializer(solventRsw_);
         serializer(micp_);
         serializer(mixControls_);
     }
@@ -360,6 +366,7 @@ protected:
     std::vector<Scalar> minOilPressure_;
     std::vector<Scalar> overburdenPressure_;
     std::vector<Scalar> solventSaturation_;
+    std::vector<Scalar> solventRsw_;
     MICPSolutionContainer<Scalar> micp_;
 
     EclMixingRateControls<FluidSystem, Scalar> mixControls_;

ebos/eclgenericproblem_impl.hh

@@ -106,6 +106,7 @@ serializationTestObject(const EclipseState& eclState,
     result.minOilPressure_ = {7.0, 8.0, 9.0, 10.0};
     result.overburdenPressure_ = {11.0};
     result.solventSaturation_ = {15.0};
+    result.solventRsw_ = {18.0};
     result.polymer_ = PolymerSolutionContainer<Scalar>::serializationTestObject();
     result.micp_ = MICPSolutionContainer<Scalar>::serializationTestObject();
     result.mixControls_ = EclMixingRateControls<FluidSystem,Scalar>::serializationTestObject(schedule);
@@ -529,6 +530,11 @@ readBlackoilExtentionsInitialConditions_(std::size_t numDof,
             solventSaturation_ = eclState_.fieldProps().get_double("SSOL");
         else
             solventSaturation_.resize(numDof, 0.0);
+
+        //if (eclState_.fieldProps().has_double("SSOL"))
+        //    solventRsw_ = eclState_.fieldProps().get_double("SSOL");
+        //else
+            solventRsw_.resize(numDof, 0.0);
     }
 
     if (enablePolymer) {
@@ -617,6 +623,17 @@ solventSaturation(unsigned elemIdx) const
     return solventSaturation_[elemIdx];
 }
 
+
+template<class GridView, class FluidSystem, class Scalar>
+Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
+solventRsw(unsigned elemIdx) const
+{
+    if (solventRsw_.empty())
+        return 0;
+
+    return solventRsw_[elemIdx];
+}
+
 template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 drsdtcon(unsigned elemIdx, int episodeIdx) const
@@ -761,6 +778,7 @@ operator==(const EclGenericProblem& rhs) const
            this->minOilPressure_ == rhs.minOilPressure_ &&
            this->overburdenPressure_ == rhs.overburdenPressure_ &&
            this->solventSaturation_ == rhs.solventSaturation_ &&
+           this->solventRsw_ == rhs.solventRsw_ &&
            this->polymer_ == rhs.polymer_ &&
            this->micp_ == rhs.micp_ &&
            this->mixControls_ == rhs.mixControls_;

ebos/ecloutputblackoilmodule.hh

@@ -416,6 +416,10 @@ public:
                 this->sSol_[globalDofIdx] = intQuants.solventSaturation().value();
             }
 
+            if (!this->rswSol_.empty()) {
+                this->rswSol_[globalDofIdx] = intQuants.rsSolw().value();
+            }
+
             if (!this->cPolymer_.empty()) {
                 this->cPolymer_[globalDofIdx] = intQuants.polymerConcentration().value();
             }

ebos/eclproblem.hh

@@ -1294,8 +1294,7 @@ public:
         values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
         values.assignNaive(initialFluidStates_[globalDofIdx]);
 
-        if constexpr (enableSolvent)
+        SolventModule::assignPrimaryVars(values, this->solventSaturation_[globalDofIdx], this->solventRsw_[globalDofIdx]);
-            values[Indices::solventSaturationIdx] = this->solventSaturation_[globalDofIdx];
 
         if constexpr (enablePolymer)
             values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx];
@@ -2141,8 +2140,10 @@ protected:
 
         std::size_t numElems = this->model().numGridDof();
         initialFluidStates_.resize(numElems);
-        if constexpr (enableSolvent)
+        if constexpr (enableSolvent) {
             this->solventSaturation_.resize(numElems, 0.0);
+            this->solventRsw_.resize(numElems, 0.0);
+        }
 
         if constexpr (enablePolymer)
             this->polymer_.concentration.resize(numElems, 0.0);
@@ -2178,8 +2179,10 @@ protected:
 
                 processRestartSaturations_(elemFluidState, ssol);
 
-                if constexpr (enableSolvent)
+                if constexpr (enableSolvent) {
                     this->solventSaturation_[elemIdx] = ssol;
+                    this->solventRsw_[elemIdx] = eclWriter_->eclOutputModule().getSolventRsw(elemIdx);
+                }
             }
 
             this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());

opm/simulators/flow/priVarsPacking.hpp

@@ -33,7 +33,8 @@ namespace Opm {
             std::size_t m2 = static_cast<std::size_t>(privar.primaryVarsMeaningPressure());
             std::size_t m3 = static_cast<std::size_t>(privar.primaryVarsMeaningGas());
             std::size_t m4 = static_cast<std::size_t>(privar.primaryVarsMeaningBrine());
-            return m1 + (m2 << fbits*1) + (m3 << fbits*2) + (m4 << fbits*3);
+            std::size_t m5 = static_cast<std::size_t>(privar.primaryVarsMeaningSolvent());
+            return m1 + (m2 << fbits*1) + (m3 << fbits*2) + (m4 << fbits*3) + (m5 << fbits*4);
         }
 
         template <class PV>
@@ -43,10 +44,12 @@ namespace Opm {
             std::size_t m2 = (meanings >> fbits*1) & filter;
             std::size_t m3 = (meanings >> fbits*2) & filter;
             std::size_t m4 = (meanings >> fbits*3) & filter;
+            std::size_t m5 = (meanings >> fbits*4) & filter;
             privar.setPrimaryVarsMeaningWater(typename PV::WaterMeaning(m1));
             privar.setPrimaryVarsMeaningPressure(typename PV::PressureMeaning(m2));
             privar.setPrimaryVarsMeaningGas(typename PV::GasMeaning(m3));
             privar.setPrimaryVarsMeaningBrine(typename PV::BrineMeaning(m4));
+            privar.setPrimaryVarsMeaningSolvent(typename PV::SolventMeaning(m5));
         }
     } // namespace PVMeanings
 } // namespace Opm

tests/test_privarspacking.cpp

@@ -56,6 +56,12 @@ public:
         Disabled, // The primary variable is not used
     };
 
+    enum class SolventMeaning {
+        Ss, // solvent saturation
+        Rsolw, // dissolved solvent in water
+        Disabled, // The primary variable is not used
+    };
+
     WaterMeaning primaryVarsMeaningWater() const
     { return primaryVarsMeaningWater_; }
     void setPrimaryVarsMeaningWater(WaterMeaning newMeaning)
@@ -76,18 +82,25 @@ public:
     void setPrimaryVarsMeaningBrine(BrineMeaning newMeaning)
     { primaryVarsMeaningBrine_ = newMeaning; }
 
+    SolventMeaning primaryVarsMeaningSolvent() const
+    { return primaryVarsMeaningSolvent_; }
+    void setPrimaryVarsMeaningSolvent(SolventMeaning newMeaning)
+    { primaryVarsMeaningSolvent_ = newMeaning; }
+
     bool operator==(const PriVarMeaning& other) const
     {
         return primaryVarsMeaningWater_ == other.primaryVarsMeaningWater_
             && primaryVarsMeaningPressure_ == other.primaryVarsMeaningPressure_
             && primaryVarsMeaningGas_ == other.primaryVarsMeaningGas_
-            && primaryVarsMeaningBrine_ == other.primaryVarsMeaningBrine_;
+            && primaryVarsMeaningBrine_ == other.primaryVarsMeaningBrine_
+            && primaryVarsMeaningSolvent_ == other.primaryVarsMeaningSolvent_;
     }
 private:
     WaterMeaning primaryVarsMeaningWater_ = WaterMeaning::Disabled;
     PressureMeaning primaryVarsMeaningPressure_ = PressureMeaning::Pw;
     GasMeaning primaryVarsMeaningGas_ = GasMeaning::Disabled;
     BrineMeaning primaryVarsMeaningBrine_ = BrineMeaning::Disabled;
+    SolventMeaning primaryVarsMeaningSolvent_ = SolventMeaning::Disabled;
 };
 
 
@@ -109,6 +122,7 @@ BOOST_AUTO_TEST_CASE(meanings)
         pv1.setPrimaryVarsMeaningWater(PriVarMeaning::WaterMeaning::Rvw);
         pv1.setPrimaryVarsMeaningGas(PriVarMeaning::GasMeaning::Disabled);
         pv1.setPrimaryVarsMeaningBrine(PriVarMeaning::BrineMeaning::Cs);
+        pv1.setPrimaryVarsMeaningSolvent(PriVarMeaning::SolventMeaning::Disabled);
         std::size_t p = Opm::PVUtil::pack(pv1);
         Opm::PVUtil::unPack(pv2, p);
         BOOST_CHECK(pv1 == pv2);