Testing with different producer switching logic including check for control feasibility

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1385,9 +1385,9 @@ updateAndCommunicateGroupData(const int reportStepIdx,
             calcResvCoeff(fipnum, pvtreg, this->groupState().production_rates(group.name()), resv_coeff);
             const Scalar efficiencyFactor = well->wellEcl().getEfficiencyFactor() *
                                     ws.efficiency_scaling_factor;
-            // Translate injector type from control to Phase.
-            Scalar group_target = std::numeric_limits<Scalar>::max();
             if (well->isProducer()) {
+                std::pair<WellProducerCMode, Scalar> group_target;
+                group_target.second = std::numeric_limits<Scalar>::max();
                 group_target = wg_helper.getWellGroupTargetProducer(
                     well->name(),
                     well->wellEcl().groupName(),
@@ -1397,6 +1397,9 @@ updateAndCommunicateGroupData(const int reportStepIdx,
                     resv_coeff,
                     deferred_logger
                 );
+                auto& ws_update = this->wellState().well(well->indexOfWell());
+                ws_update.production_cmode_group_translated = group_target.first;
+                ws_update.group_target = group_target.second;
             } else {
                 const auto& well_controls = well->wellEcl().injectionControls(summaryState_);
                 auto injectorType = well_controls.injector_type;
@@ -1420,6 +1423,7 @@ updateAndCommunicateGroupData(const int reportStepIdx,
                 default:
                     throw std::logic_error("MULTI-phase injection is not supported, but was requested for well " + well->name());
                 }
+                std::pair<WellInjectorCMode, Scalar> group_target;
                 group_target = wg_helper.getWellGroupTargetInjector(
                     well->name(),
                     well->wellEcl().groupName(),
@@ -1430,9 +1434,10 @@ updateAndCommunicateGroupData(const int reportStepIdx,
                     resv_coeff,
                     deferred_logger
                 );
+                auto& ws_update = this->wellState().well(well->indexOfWell());
+                ws_update.injection_cmode_group_translated = group_target.first;
+                ws_update.group_target = group_target.second;
             }
-            auto& ws_update = this->wellState().well(well->indexOfWell());
-            ws_update.group_target = group_target;
         }
     }
 }

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -681,6 +681,53 @@ getWQTotal() const
     return evaluation_[0][WQTotal];
 }
 
+template<class FluidSystem, class Indices>
+void MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
+fetchWellSurfaceFractions(std::vector<Scalar>& surface_fractions) const
+{
+    surface_fractions.resize(well_.numPhases(), 0.0);
+    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+        const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
+        surface_fractions[oil_pos] = 1.0;
+        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+            surface_fractions[water_pos] = value_[0][WFrac];
+            surface_fractions[oil_pos] -= surface_fractions[water_pos];
+        }
+
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            surface_fractions[gas_pos] = value_[0][GFrac];
+            surface_fractions[oil_pos] -= surface_fractions[gas_pos];
+        }
+    }
+    else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+        const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+        surface_fractions[water_pos] = 1.0;
+
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            surface_fractions[gas_pos] = value_[0][GFrac];
+            surface_fractions[water_pos] -= surface_fractions[gas_pos];
+        }
+    }
+    else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+        const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+        surface_fractions[gas_pos] = 1.0;
+    }
+    for (int p = 0; p < well_.numPhases(); ++p) {
+        const Scalar scale = well_.scalingFactor(p);
+        // for injection wells, there should only one non-zero scaling factor
+        if (scale > 0.) {
+            surface_fractions[p] /= scale;
+        } else {
+            // this should only happen to injection wells
+            surface_fractions[p] = 0.;
+        }
+        surface_fractions[p] = std::max(surface_fractions[p], Scalar{0.0});
+    }
+}
+
 template<class FluidSystem, class Indices>
 void
 MultisegmentWellPrimaryVariables<FluidSystem,Indices>::

opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp

@@ -143,6 +143,9 @@ class MultisegmentWellPrimaryVariables
     void setValue(const int idx, const std::array<Scalar, numWellEq>& val)
     { value_[idx] = val; }
 
+    //! \brief Fetch surface fractions for the well.
+    void fetchWellSurfaceFractions(std::vector<Scalar> &surface_fractions) const;
+
     //! output the segments with pressure close to lower pressure limit for debugging purpose
     void outputLowLimitPressureSegments(DeferredLogger& deferred_logger) const;
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -44,6 +44,7 @@
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/simulators/wells/ParallelWellInfo.hpp>
+#include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
 
 #include <algorithm>
 #include <cstddef>
@@ -1638,8 +1639,6 @@ namespace Opm
                 return false;
         }
 
-        updatePrimaryVariables(simulator, well_state, deferred_logger);
-
         std::vector<std::vector<Scalar> > residual_history;
         std::vector<Scalar> measure_history;
         int it = 0;
@@ -1664,22 +1663,35 @@ namespace Opm
         const auto operability_orig = this->operability_status_;
         auto well_status_cur = well_status_orig;
         // don't allow opening wells that has a stopped well status
-        const bool allow_open = well_state.well(this->index_of_well_).status == WellStatus::OPEN;
+        auto& ws = well_state.well(this->index_of_well_); 
+        const bool allow_open = ws.status == WellStatus::OPEN;
         // don't allow switcing for wells under zero rate target or requested fixed status and control
         const bool allow_switching = !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
                                      (!fixed_control || !fixed_status) && allow_open;
         bool final_check = false;
         // well needs to be set operable or else solving/updating of re-opened wells is skipped
         this->operability_status_.resetOperability();
         this->operability_status_.solvable = true;
+        updatePrimaryVariables(simulator, well_state, deferred_logger);
+        // keep track of surface rate fractions for checking control feasibility
+        std::vector<Scalar> surface_fractions;
+        bool feasible_constraint = true;
 
         for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
             ++its_since_last_switch;
-            if (allow_switching && its_since_last_switch >= min_its_after_switch && status_switch_count < max_status_switch){
+            if (this->isProducer()) {
+                this->primary_variables_.fetchWellSurfaceFractions(surface_fractions);
+                feasible_constraint = this->isFeasibleProductionConstraint(ws, surface_fractions, ws.production_cmode);
+            }
+            const bool check_controls = allow_switching &&
+                                        its_since_last_switch >= min_its_after_switch &&
+                                        status_switch_count < max_status_switch;
+            // if constraint is not feasible, we force a switch to prevent singularity
+            if (check_controls || !feasible_constraint) {
                 const Scalar wqTotal = this->primary_variables_.getWQTotal().value();
                 bool changed = this->updateWellControlAndStatusLocalIteration(
-                    simulator, wgHelper, inj_controls, prod_controls, wqTotal,
-                    well_state, deferred_logger, fixed_control, fixed_status
+                    simulator, wgHelper, inj_controls, prod_controls, wqTotal, well_state,
+                    surface_fractions, deferred_logger, fixed_control, fixed_status
                 );
                 if (changed) {
                     its_since_last_switch = 0;

opm/simulators/wells/SingleWellState.cpp

@@ -53,6 +53,7 @@ SingleWellState(const std::string& name_,
     , prev_surface_rates(pu.numActivePhases())
     , perf_data(perf_input.size(), pressure_first_connection, !is_producer, pu.numActivePhases())
     , trivial_group_target(false)
+    , prevent_group_control(false)
 {
     for (std::size_t perf = 0; perf < perf_input.size(); perf++) {
         this->perf_data.cell_index[perf] = perf_input[perf].cell_index;
@@ -410,7 +411,10 @@ bool SingleWellState<Scalar, IndexTraits>::operator==(const SingleWellState& rhs
            this->production_cmode == rhs.production_cmode &&
            this->alq_state == rhs.alq_state &&
            this->primaryvar == rhs.primaryvar &&
-           this->group_target == rhs.group_target;
+           this->group_target == rhs.group_target &&
+           this->injection_cmode_group_translated == rhs.injection_cmode_group_translated &&
+           this->production_cmode_group_translated == rhs.production_cmode_group_translated &&
+           this->prevent_group_control == rhs.prevent_group_control;
 }
 
 template class SingleWellState<double, BlackOilDefaultFluidSystemIndices>;

opm/simulators/wells/SingleWellState.hpp

@@ -85,6 +85,9 @@ class SingleWellState {
         serializer(primaryvar);
         serializer(alq_state);
         serializer(group_target);
+        serializer(injection_cmode_group_translated);
+        serializer(production_cmode_group_translated);
+        serializer(prevent_group_control);
     }
 
     bool operator==(const SingleWellState&) const;
@@ -125,6 +128,9 @@ class SingleWellState {
     PerfData<Scalar> perf_data;
     bool trivial_group_target;
     std::optional<Scalar> group_target;
+    std::optional<WellInjectorCMode> injection_cmode_group_translated{WellInjectorCMode::CMODE_UNDEFINED};
+    std::optional<WellProducerCMode> production_cmode_group_translated{WellProducerCMode::CMODE_UNDEFINED};
+    bool prevent_group_control;
     SegmentState<Scalar> segments;
     Events events;
     WellInjectorCMode injection_cmode{WellInjectorCMode::CMODE_UNDEFINED};

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -760,6 +760,53 @@ checkFinite(DeferredLogger& deferred_logger) const
     }
 }
 
+template<class FluidSystem, class Indices>
+void StandardWellPrimaryVariables<FluidSystem,Indices>::
+fetchWellSurfaceFractions(std::vector<Scalar>& surface_fractions) const
+{
+    surface_fractions.resize(well_.numPhases(), 0.0);
+    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+        const int oil_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::oilPhaseIdx);
+        surface_fractions[oil_pos] = 1.0;
+        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+            surface_fractions[water_pos] = value_[WFrac];
+            surface_fractions[oil_pos] -= surface_fractions[water_pos];
+        }
+
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            surface_fractions[gas_pos] = value_[GFrac];
+            surface_fractions[oil_pos] -= surface_fractions[gas_pos];
+        }
+    }
+    else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+        const int water_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::waterPhaseIdx);
+        surface_fractions[water_pos] = 1.0;
+
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+            surface_fractions[gas_pos] = value_[GFrac];
+            surface_fractions[water_pos] -= surface_fractions[gas_pos];
+        }
+    }
+    else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+        const int gas_pos = FluidSystem::canonicalToActivePhaseIdx(FluidSystem::gasPhaseIdx);
+        surface_fractions[gas_pos] = 1.0;
+    }
+    for (int p = 0; p < well_.numPhases(); ++p) {
+        const Scalar scale = well_.scalingFactor(p);
+        // for injection wells, there should be only one non-zero scaling factor
+        if (scale > 0.) {
+            surface_fractions[p] /= scale;
+        } else {
+            // this should only happen to injection wells
+            surface_fractions[p] = 0.;
+        }
+        surface_fractions[p] = std::max(surface_fractions[p], Scalar{0.0});
+    }
+}
+
 #include <opm/simulators/utils/InstantiationIndicesMacros.hpp>
 
 INSTANTIATE_TYPE_INDICES(StandardWellPrimaryVariables, double)

opm/simulators/wells/StandardWellPrimaryVariables.hpp

@@ -148,6 +148,9 @@ class StandardWellPrimaryVariables {
     void setValue(const int idx, const Scalar val)
     { value_[idx] = val; }
 
+    //! \brief Fetch well surface fraction values.
+    void fetchWellSurfaceFractions(std::vector<Scalar>& surface_fractions) const;
+
 private:
     //! \brief Initialize evaluations from values.
     void setEvaluationsFromValues();

opm/simulators/wells/StandardWell_impl.hpp

@@ -37,6 +37,7 @@
 #include <opm/simulators/wells/VFPHelpers.hpp>
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
 #include <opm/simulators/wells/WellConvergence.hpp>
+#include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
 
 #include <algorithm>
 #include <cstddef>
@@ -2444,7 +2445,6 @@ namespace Opm
                                const bool fixed_status /*false*/)
     {
         const auto& group_state = wgHelper.groupState();
-        updatePrimaryVariables(simulator, well_state, deferred_logger);
 
         const int max_iter = this->param_.max_inner_iter_wells_;
         int it = 0;
@@ -2467,7 +2467,8 @@ namespace Opm
         auto well_status_cur = well_status_orig;
         int status_switch_count = 0;
         // don't allow opening wells that has a stopped well status
-        const bool allow_open = well_state.well(this->index_of_well_).status == WellStatus::OPEN;
+        auto& ws = well_state.well(this->index_of_well_);
+        const bool allow_open = ws.status == WellStatus::OPEN;
         // don't allow switcing for wells under zero rate target or requested fixed status and control
         const bool allow_switching =
             !this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) &&
@@ -2478,13 +2479,25 @@ namespace Opm
         // well needs to be set operable or else solving/updating of re-opened wells is skipped
         this->operability_status_.resetOperability();
         this->operability_status_.solvable = true;
+        updatePrimaryVariables(simulator, well_state, deferred_logger);
+        // keep track of surface rate fractions for checking control feasibility
+        std::vector<Scalar> surface_fractions;
+        bool feasible_constraint = true;
         do {
             its_since_last_switch++;
-            if (allow_switching && its_since_last_switch >= min_its_after_switch && status_switch_count < max_status_switch){
+            if (this->isProducer()) {
+                this->primary_variables_.fetchWellSurfaceFractions(surface_fractions);
+                feasible_constraint = this->isFeasibleProductionConstraint(ws, surface_fractions, ws.production_cmode);
+            }
+            const bool check_controls = allow_switching &&
+                                        its_since_last_switch >= min_its_after_switch &&
+                                        status_switch_count < max_status_switch;
+            // if constraint is not feasible, we force a switch to prevent singularity
+            if (check_controls || !feasible_constraint){
                 const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value();
                 changed = this->updateWellControlAndStatusLocalIteration(
-                    simulator, wgHelper, inj_controls, prod_controls, wqTotal,
-                    well_state, deferred_logger, fixed_control, fixed_status
+                    simulator, wgHelper, inj_controls, prod_controls, wqTotal, well_state,
+                    surface_fractions, deferred_logger, fixed_control, fixed_status
                 );
                 if (changed){
                     its_since_last_switch = 0;