Merge pull request #4549 from totto82/resv_undersat

fix RESV for undersaturated wells

opm/simulators/wells/BlackoilWellModel.hpp

@@ -444,6 +444,7 @@ namespace Opm {
 
             void calcRates(const int fipnum,
                            const int pvtreg,
+                           const std::vector<double>& production_rates,
                            std::vector<double>& resv_coeff) override;
 
             void calcInjRates(const int fipnum,

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -549,7 +549,7 @@ checkGroupHigherConstraints(const Group& group,
             rates[phasePos] = -comm_.sum(local_current_rate);
         }
         std::vector<double> resv_coeff(phase_usage_.num_phases, 0.0);
-        calcRates(fipnum, pvtreg, resv_coeff);
+        calcRates(fipnum, pvtreg, this->groupState().production_rates(group.name()), resv_coeff);
         // Check higher up only if under individual (not FLD) control.
         const Group::ProductionCMode& currentControl = this->groupState().production_control(group.name());
         if (currentControl != Group::ProductionCMode::FLD && group.productionGroupControlAvailable()) {

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -312,6 +312,7 @@ protected:
                      double wsolvent);
     virtual void calcRates(const int fipnum,
                            const int pvtreg,
+                           const std::vector<double>& production_rates,
                            std::vector<double>& resv_coeff) = 0;
     virtual void calcInjRates(const int fipnum,
                            const int pvtreg,

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1859,9 +1859,10 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     calcRates(const int fipnum,
               const int pvtreg,
+              const std::vector<double>& production_rates,
               std::vector<double>& resv_coeff)
     {
-        rateConverter_->calcCoeff(fipnum, pvtreg, resv_coeff);
+        rateConverter_->calcCoeff(fipnum, pvtreg, production_rates, resv_coeff);
     }
 
     template<typename TypeTag>

opm/simulators/wells/RateConverter.cpp

@@ -128,11 +128,43 @@ template <class Coeff>
 void SurfaceToReservoirVoidage<FluidSystem,Region>::
 calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
 {
-    const auto& pu = phaseUsage_;
     const auto& ra = attr_.attributes(r);
-    const double p = ra.pressure;
+    calcCoeff(pvtRegionIdx, ra.pressure, ra.rs, ra.rv, ra.rsw, ra.rvw, ra.temperature, ra.saltConcentration, coeff);
-    const double T = ra.temperature;
+}
-    const double saltConcentration = ra.saltConcentration;
+
+template <class FluidSystem, class Region>
+template <class Coeff, class Rates>
+void SurfaceToReservoirVoidage<FluidSystem,Region>::
+calcCoeff(const RegionId r, const int pvtRegionIdx, const Rates& surface_rates, Coeff& coeff) const
+{
+    const auto& ra = attr_.attributes(r);
+    const auto& pu = phaseUsage_;
+    const int   iw = RegionAttributeHelpers::PhasePos::water(pu);
+    const int   io = RegionAttributeHelpers::PhasePos::oil  (pu);
+    const int   ig = RegionAttributeHelpers::PhasePos::gas  (pu);
+    const auto [Rs, Rv] =
+        dissolvedVaporisedRatio(io, ig, ra.rs, ra.rv, surface_rates);
+
+    const auto [Rsw, Rvw] =
+        dissolvedVaporisedRatio(iw, ig, ra.rsw, ra.rvw, surface_rates);
+
+    calcCoeff(pvtRegionIdx, ra.pressure, Rs, Rv, Rsw, Rvw, ra.temperature, ra.saltConcentration, coeff);
+}
+
+template <class FluidSystem, class Region>
+template <class Coeff>
+void SurfaceToReservoirVoidage<FluidSystem,Region>::
+calcCoeff(  const int pvtRegionIdx,
+            const double        p,
+            const double        Rs,
+            const double        Rv,
+            const double        Rsw,
+            const double        Rvw,
+            const double        T,
+            const double saltConcentration,
+            Coeff& coeff) const
+{
+    const auto& pu = phaseUsage_;
 
     const int   iw = RegionAttributeHelpers::PhasePos::water(pu);
     const int   io = RegionAttributeHelpers::PhasePos::oil  (pu);
@@ -140,9 +172,6 @@ calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
 
     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);
 
@@ -156,21 +185,19 @@ calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
         coeff[iw] += 1.0 / den;
 
         if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
-            coeff[ig] -= ra.rvw / den;
+            coeff[ig] -= Rvw / den;
         }
     }
 
-    // Actual Rs and Rv:
-    double Rs = ra.rs;
-    double Rv = ra.rv;
-
     // 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);
+        std::string msg = "We currently support either gas in water or gas in oil, not both."
+        "i.e. detR = " + std::to_string(detR) + " and detRw " + std::to_string(detRw) +
+        "can not both be different from 1";
         throw std::range_error(msg);
     }
 
@@ -183,7 +210,7 @@ calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
         coeff[io] += 1.0 / den;
 
         if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
-            coeff[ig] -= ra.rv / den;
+            coeff[ig] -= Rv / den;
         }
     }
 
@@ -194,18 +221,19 @@ calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
             const double denw = bg * detRw;
             coeff[ig] += 1.0 / denw;
             if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
-               coeff[iw] -= ra.rsw / denw;
+               coeff[iw] -= Rsw / denw;
             }
         } else {
             const double den = bg * detR;
             coeff[ig] += 1.0 / den;
             if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
-                coeff[io] -= ra.rs / den;
+                coeff[io] -= Rs / den;
             }
         }
     }
 }
 
+
 template <class FluidSystem, class Region>
 template <typename SurfaceRates, typename VoidageRates>
 void SurfaceToReservoirVoidage<FluidSystem,Region>::
@@ -340,6 +368,8 @@ template void SurfaceToReservoirVoidage<FS,std::vector<int>>::
     calcInjCoeff<std::vector<double>>(const int, const int, std::vector<double>&) const;
 template void SurfaceToReservoirVoidage<FS,std::vector<int>>::
     calcCoeff<std::vector<double>>(const int, const int, std::vector<double>&) const;
+template void SurfaceToReservoirVoidage<FS,std::vector<int>>::
+    calcCoeff<std::vector<double>, std::vector<double>>(const int, const int, const std::vector<double>&, std::vector<double>&) const;
 template void SurfaceToReservoirVoidage<FS,std::vector<int>>::
     calcReservoirVoidageRates<std::vector<double>,std::vector<double>>(const int,
                                                                        const double,

opm/simulators/wells/RateConverter.hpp

@@ -243,6 +243,22 @@ namespace Opm {
             void
             calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const;
 
+            template <class Coeff , class Rates>
+            void
+            calcCoeff(const RegionId r, const int pvtRegionIdx, const Rates& surface_rates, Coeff& coeff) const;
+
+            template <class Coeff>
+            void
+            calcCoeff(  const int pvtRegionIdx,
+                        const double        p,
+                        const double        rs,
+                        const double        rv,
+                        const double        rsw,
+                        const double        rvw,
+                        const double        T,
+                        const double saltConcentration,
+                        Coeff& coeff) const;
+
             template <class Coeff>
             void
             calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const;

opm/simulators/wells/WellAssemble.cpp

@@ -103,7 +103,7 @@ assembleControlEqProd(const WellState& well_state,
         auto total_rate = rates[0]; // To get the correct type only.
         total_rate = 0.0;
         std::vector<double> convert_coeff(well_.numPhases(), 1.0);
-        well_.rateConverter().calcCoeff(/*fipreg*/ 0, well_.pvtRegionIdx(), convert_coeff);
+        well_.rateConverter().calcCoeff(/*fipreg*/ 0, well_.pvtRegionIdx(), well_state.well(well_.indexOfWell()).surface_rates, convert_coeff);
         for (int phase = 0; phase < 3; ++phase) {
             if (pu.phase_used[phase]) {
                 const int pos = pu.phase_pos[phase];
@@ -153,9 +153,13 @@ assembleControlEqProd(const WellState& well_state,
                 active_rates[pu.phase_pos[canonical_phase]] = rates[canonical_phase];
             }
         }
-        auto rCoeff = [this](const RegionId id, const int region, std::vector<double>& coeff)
+        auto rCoeff = [this, &group_state](const RegionId id, const int region, const std::optional<std::string>& prod_gname, std::vector<double>& coeff)
         {
-            well_.rateConverter().calcCoeff(id, region, coeff);
+            if (prod_gname)
+                well_.rateConverter().calcCoeff(id, region, group_state.production_rates(*prod_gname), coeff);
+            else
+                well_.rateConverter().calcCoeff(id, region, coeff);
+
         };
         WellGroupControls(well_).getGroupProductionControl(group, well_state,
                                                              group_state,
@@ -243,9 +247,13 @@ assembleControlEqInj(const WellState& well_state,
     case Well::InjectorCMode::GRUP: {
         assert(well_.wellEcl().isAvailableForGroupControl());
         const auto& group = schedule.getGroup(well_.wellEcl().groupName(), well_.currentStep());
-        auto rCoeff = [this](const RegionId id, const int region, std::vector<double>& coeff)
+        auto rCoeff = [this, &group_state](const RegionId id, const int region, const std::optional<std::string>& prod_gname, std::vector<double>& coeff)
         {
-            well_.rateConverter().calcInjCoeff(id, region, coeff);
+            if(prod_gname) {
+                well_.rateConverter().calcCoeff(id, region, group_state.production_rates(*prod_gname), coeff);
+            } else {
+                well_.rateConverter().calcInjCoeff(id, region, coeff);
+            }
         };
         WellGroupControls(well_).getGroupInjectionControl(group,
                                                             well_state,

opm/simulators/wells/WellGroupConstraints.cpp

@@ -70,7 +70,7 @@ checkGroupConstraintsInj(const Group& group,
 
     // Make conversion factors for RESV <-> surface rates.
     std::vector<double> resv_coeff(well_.phaseUsage().num_phases, 1.0);
-    rateConverter(0, well_.pvtRegionIdx(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
+    rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
 
     const auto& ws = well_state.well(well_.indexOfWell());
     // Call check for the well's injection phase.
@@ -104,7 +104,7 @@ checkGroupConstraintsProd(const Group& group,
 {
     // Make conversion factors for RESV <-> surface rates.
     std::vector<double> resv_coeff(well_.phaseUsage().num_phases, 1.0);
-    rateConverter(0, well_.pvtRegionIdx(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
+    rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
 
     const auto& ws = well_state.well(well_.indexOfWell());
     return WellGroupHelpers::checkGroupConstraintsProd(well_.name(),

opm/simulators/wells/WellGroupConstraints.hpp

@@ -27,6 +27,8 @@
 #include <functional>
 #include <utility>
 #include <vector>
+#include <string>
+#include <optional>
 
 namespace Opm
 {
@@ -47,7 +49,7 @@ public:
     //! \brief Constructor sets reference to well.
     WellGroupConstraints(const WellInterfaceGeneric& well) : well_(well) {}
 
-    using RateConvFunc = std::function<void(const RegionId, const int, std::vector<double>&)>;
+    using RateConvFunc = std::function<void(const RegionId, const int, const std::optional<std::string>&, std::vector<double>&)>;
 
     bool checkGroupConstraints(WellState& well_state,
                                const GroupState& group_state,

opm/simulators/wells/WellGroupControls.cpp

@@ -127,7 +127,7 @@ getGroupInjectionControl(const Group& group,
 
     // Make conversion factors for RESV <-> surface rates.
     std::vector<double> resv_coeff(pu.num_phases, 1.0);
-    rateConverter(0, well_.pvtRegionIdx(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
+    rateConverter(0, well_.pvtRegionIdx(), std::nullopt, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
 
     double sales_target = 0;
     if (schedule[well_.currentStep()].gconsale().has(group.name())) {
@@ -251,7 +251,7 @@ getGroupInjectionTargetRate(const Group& group,
 
     // Make conversion factors for RESV <-> surface rates.
     std::vector<double> resv_coeff(pu.num_phases, 1.0);
-    rateConverter(0, well_.pvtRegionIdx(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
+    rateConverter(0, well_.pvtRegionIdx(), std::nullopt, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
 
     double sales_target = 0;
     if (schedule[well_.currentStep()].gconsale().has(group.name())) {
@@ -353,7 +353,7 @@ void WellGroupControls::getGroupProductionControl(const Group& group,
 
     // Make conversion factors for RESV <-> surface rates.
     std::vector<double> resv_coeff(well_.phaseUsage().num_phases, 1.0);
-    rateConverter(0, well_.pvtRegionIdx(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
+    rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
 
     // gconsale may adjust the grat target.
     // the adjusted rates is send to the targetCalculator
@@ -442,7 +442,7 @@ getGroupProductionTargetRate(const Group& group,
 
     // Make conversion factors for RESV <-> surface rates.
     std::vector<double> resv_coeff(well_.phaseUsage().num_phases, 1.0);
-    rateConverter(0, well_.pvtRegionIdx(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
+    rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
 
     // gconsale may adjust the grat target.
     // the adjusted rates is send to the targetCalculator

opm/simulators/wells/WellGroupControls.hpp

@@ -24,6 +24,7 @@
 #ifndef OPM_WELL_GROUP_CONTROLS_HEADER_INCLUDED
 #define OPM_WELL_GROUP_CONTROLS_HEADER_INCLUDED
 
+#include <string>
 #include <functional>
 #include <optional>
 #include <vector>
@@ -47,7 +48,7 @@ public:
     //! \brief Constructor sets reference to well.
     WellGroupControls(const WellInterfaceGeneric& well) : well_(well) {}
 
-    using RateConvFunc = std::function<void(const RegionId, const int, std::vector<double>&)>;
+    using RateConvFunc = std::function<void(const RegionId, const int, const std::optional<std::string>&, std::vector<double>&)>;
 
     template<class EvalWell>
     void getGroupInjectionControl(const Group& group,

opm/simulators/wells/WellInterfaceFluidSystem.cpp

@@ -131,9 +131,13 @@ checkGroupConstraints(WellState& well_state,
                       const SummaryState& summaryState,
                       DeferredLogger& deferred_logger) const
 {
-    auto rCoeff = [this](const int id, const int region, std::vector<double>& coeff)
+    auto rCoeff = [this, &group_state](const RegionId id, const int region, const std::optional<std::string>& prod_gname, std::vector<double>& coeff)
     {
-        this->rateConverter().calcCoeff(id, region, coeff);
+        if (prod_gname)
+            this->rateConverter().calcCoeff(id, region, group_state.production_rates(*prod_gname), coeff);
+        else
+            this->rateConverter().calcInjCoeff(id, region, coeff);
+
     };
 
     return WellGroupConstraints(*this).checkGroupConstraints(well_state, group_state,
@@ -187,9 +191,13 @@ getGroupInjectionTargetRate(const Group& group,
                             double efficiencyFactor,
                             DeferredLogger& deferred_logger) const
 {
-    auto rCoeff = [this](const int id, const int region, std::vector<double>& coeff)
+    auto rCoeff = [this, &group_state](const RegionId id, const int region, const std::optional<std::string>& prod_gname, std::vector<double>& coeff)
     {
-        this->rateConverter().calcCoeff(id, region, coeff);
+        if (prod_gname)
+            this->rateConverter().calcCoeff(id, region, group_state.production_rates(*prod_gname), coeff);
+        else
+            this->rateConverter().calcInjCoeff(id, region, coeff);
+
     };
 
     return WellGroupControls(*this).getGroupInjectionTargetRate(group, well_state,
@@ -210,9 +218,13 @@ getGroupProductionTargetRate(const Group& group,
                           double efficiencyFactor,
                           DeferredLogger& deferred_logger) const
 {
-    auto rCoeff = [this](const int id, const int region, std::vector<double>& coeff)
+    auto rCoeff = [this, &group_state](const RegionId id, const int region, const std::optional<std::string>& prod_gname, std::vector<double>& coeff)
     {
-        this->rateConverter().calcCoeff(id, region, coeff);
+        if (prod_gname)
+            this->rateConverter().calcCoeff(id, region, group_state.production_rates(*prod_gname), coeff);
+        else
+            this->rateConverter().calcInjCoeff(id, region, coeff);
+
     };
 
     return WellGroupControls(*this).getGroupProductionTargetRate(group, well_state,

opm/simulators/wells/WellInterface_impl.hpp

@@ -900,7 +900,7 @@ namespace Opm
             case Well::ProducerCMode::RESV:
             {
                 std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
-                this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
+                this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, ws.surface_rates, convert_coeff);
                 double total_res_rate = 0.0;
                 for (int p = 0; p<np; ++p) {
                     total_res_rate -= ws.surface_rates[p] * convert_coeff[p];

tests/test_RestartSerialization.cpp

@@ -275,7 +275,7 @@ public:
         switched_inj_groups_ = {{{"test4", Phase::SOLVENT}, "test5"}};
     }
 
-    void calcRates(const int, const int, std::vector<double>&) override
+    void calcRates(const int, const int, const std::vector<double>&, std::vector<double>&) override
     {}
 
     void calcInjRates(const int, const int, std::vector<double>&) override