Merge pull request #4563 from GitPaean/treating_zero_injection_rate_well_as_stoppped

for the wells under zero rate target control, the THP output will be zero

opm/simulators/wells/MultisegmentWellAssemble.cpp

@@ -116,7 +116,7 @@ assembleControlEq(const WellState& well_state,
         return rates;
     };
 
-    if (well_.wellIsStopped()) {
+    if (well_.stopppedOrZeroRateTarget(summaryState, well_state)) {
         control_eq = primary_variables.getWQTotal();
     } else if (well_.isInjector() ) {
         // Find scaling factor to get injection rate,
@@ -164,9 +164,6 @@ assembleControlEq(const WellState& well_state,
                                                  bhp_from_thp,
                                                  control_eq,
                                                  deferred_logger);
-    } else if (wellhelpers::rateControlWithZeroTarget(well_state.well(well_.indexOfWell()).production_cmode, prod_controls)) {
-        // Production mode, zero target. Treat as STOP.
-        control_eq = primary_variables.getWQTotal();
     } else {
         // Find rates.
         const auto rates = getRates();

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -64,7 +64,7 @@ init()
 
 template<class FluidSystem, class Indices, class Scalar>
 void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
-update(const WellState& well_state, const bool zero_rate_target)
+update(const WellState& well_state, const bool stop_or_zero_rate_target)
 {
     static constexpr int Water = BlackoilPhases::Aqua;
     static constexpr int Gas = BlackoilPhases::Vapour;
@@ -104,7 +104,7 @@ update(const WellState& well_state, const bool zero_rate_target)
             }
         }
         value_[seg][WQTotal] = total_seg_rate;
-        if (zero_rate_target && seg == 0) {
+        if (stop_or_zero_rate_target && seg == 0) {
             value_[seg][WQTotal] = 0;
         }
         if (std::abs(total_seg_rate) > 0.) {
@@ -155,7 +155,7 @@ void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
 updateNewton(const BVectorWell& dwells,
              const double relaxation_factor,
              const double dFLimit,
-             const bool zero_rate_target,
+             const bool stop_or_zero_rate_target,
              const double max_pressure_change)
 {
     const std::vector<std::array<double, numWellEq>> old_primary_variables = value_;
@@ -198,7 +198,7 @@ updateNewton(const BVectorWell& dwells,
         }
     }
 
-    if (zero_rate_target) {
+    if (stop_or_zero_rate_target) {
         value_[0][WQTotal] = 0.;
     }
 }
@@ -207,6 +207,7 @@ template<class FluidSystem, class Indices, class Scalar>
 void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
 copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
                 const double rho,
+                const bool stop_or_zero_rate_target,
                 WellState& well_state,
                 DeferredLogger& deferred_logger) const
 {
@@ -393,7 +394,7 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
     }
 
     WellBhpThpCalculator(well_)
-        .updateThp(rho, [this]() { return well_.wellEcl().alq_value(); },
+        .updateThp(rho, stop_or_zero_rate_target, [this]() { return well_.wellEcl().alq_value(); },
                    {FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
                     FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),
                     FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)},

opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp

@@ -88,18 +88,19 @@ public:
     void init();
 
     //! \brief Copy values from well state.
-    void update(const WellState& well_state, const bool zero_rate_target);
+    void update(const WellState& well_state, const bool stop_or_zero_rate_target);
 
     //! \brief Update values from newton update vector.
     void updateNewton(const BVectorWell& dwells,
                       const double relaxation_factor,
                       const double DFLimit,
-                      const bool zero_rate_target,
+                      const bool stop_or_zero_rate_target,
                       const double max_pressure_change);
 
     //! \brief Copy values to well state.
     void copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
                          const double rho,
+                         const bool stop_or_zero_rate_target,
                          WellState& well_state,
                          DeferredLogger& deferred_logger) const;
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -152,8 +152,8 @@ namespace Opm
                            const WellState& well_state,
                            DeferredLogger& /* deferred_logger */)
     {
-        const bool zero_rate_target = this->wellUnderZeroProductionRateControl(summary_state, well_state);
+        const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
-        this->primary_variables_.update(well_state, zero_rate_target);
+        this->primary_variables_.update(well_state, stop_or_zero_rate_target);
     }
 
 
@@ -615,14 +615,14 @@ namespace Opm
 
         const double dFLimit = this->param_.dwell_fraction_max_;
         const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
-        const bool zero_rate_target = this->wellUnderZeroProductionRateControl(summary_state, well_state);
+        const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
         this->primary_variables_.updateNewton(dwells,
                                               relaxation_factor,
                                               dFLimit,
-                                              zero_rate_target,
+                                              stop_or_zero_rate_target,
                                               max_pressure_change);
 
-        this->primary_variables_.copyToWellState(*this, getRefDensity(),
+        this->primary_variables_.copyToWellState(*this, getRefDensity(), stop_or_zero_rate_target,
                                                  well_state, deferred_logger);
         Base::calculateReservoirRates(well_state.well(this->index_of_well_));
     }

opm/simulators/wells/StandardWell.hpp

@@ -365,10 +365,12 @@ namespace Opm
                                             DeferredLogger& deferred_logger) const;
 
         void updatePrimaryVariablesNewton(const BVectorWell& dwells,
-                                          const bool zero_rate_target,
+                                          const bool stop_or_zero_rate_target,
                                           DeferredLogger& deferred_logger);
 
-        void updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger) const;
+        void updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
+                                                 WellState& well_state,
+                                                 DeferredLogger& deferred_logger) const;
 
         virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
                                                     const double dt,

opm/simulators/wells/StandardWellAssemble.cpp

@@ -116,7 +116,7 @@ assembleControlEq(const WellState& well_state,
         return rates;
     };
 
-    if (well_.wellIsStopped()) {
+    if (well_.stopppedOrZeroRateTarget(summaryState, well_state)) {
         control_eq = primary_variables.eval(PrimaryVariables::WQTotal);
     } else if (well_.isInjector()) {
         // Find injection rate.
@@ -143,9 +143,6 @@ assembleControlEq(const WellState& well_state,
                                  bhp_from_thp,
                                  control_eq,
                                  deferred_logger);
-    } else if (wellhelpers::rateControlWithZeroTarget(well_state.well(well_.indexOfWell()).production_cmode, prod_controls)) {
-        // Production mode, zero target. Treat as STOP.
-        control_eq = primary_variables.eval(PrimaryVariables::WQTotal);
     } else {
              // Find rates.
         const auto rates = getRates();

opm/simulators/wells/StandardWellEval.cpp

@@ -71,13 +71,15 @@ extendEval(const Eval& in) const
 template<class FluidSystem, class Indices, class Scalar>
 void
 StandardWellEval<FluidSystem,Indices,Scalar>::
-updateWellStateFromPrimaryVariables(WellState& well_state,
+updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
+                                    WellState& well_state,
                                     DeferredLogger& deferred_logger) const
 {
     this->primary_variables_.copyToWellState(well_state, deferred_logger);
 
     WellBhpThpCalculator(baseif_).
             updateThp(connections_.rho(),
+                      stop_or_zero_rate_target,
                       [this,&well_state]() { return this->baseif_.getALQ(well_state); },
                       {FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
                        FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),

opm/simulators/wells/StandardWellEval.hpp

@@ -92,7 +92,8 @@ protected:
               const int num_cells,
               const bool has_polymermw);
 
-    void updateWellStateFromPrimaryVariables(WellState& well_state,
+    void updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
+                                             WellState& well_state,
                                              DeferredLogger& deferred_logger) const;
 
     PrimaryVariables primary_variables_; //!< Primary variables for well

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -122,7 +122,7 @@ resize(const int numWellEq)
 template<class FluidSystem, class Indices, class Scalar>
 void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
 update(const WellState& well_state,
-       const bool zero_rate_target,
+       const bool stop_or_zero_rate_target,
        DeferredLogger& deferred_logger)
 {
     static constexpr int Water = BlackoilPhases::Aqua;
@@ -160,7 +160,7 @@ update(const WellState& well_state,
         }
     } else {
             value_[WQTotal] = total_well_rate;
-            if (zero_rate_target) {
+            if (stop_or_zero_rate_target) {
                 value_[WQTotal] = 0.;
             }
     }
@@ -245,7 +245,7 @@ updatePolyMW(const WellState& well_state)
 template<class FluidSystem, class Indices, class Scalar>
 void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
 updateNewton(const BVectorWell& dwells,
-             const bool zero_rate_target,
+             const bool stop_or_zero_rate_target,
              [[maybe_unused]] const double dFLimit,
              const double dBHPLimit)
 {
@@ -281,7 +281,7 @@ updateNewton(const BVectorWell& dwells,
 
     // updating the total rates Q_t
     value_[WQTotal] = value_[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
-    if (zero_rate_target) {
+    if (stop_or_zero_rate_target) {
         value_[WQTotal] = 0.;
     }
     // TODO: here, we make sure it is zero for zero rated wells

opm/simulators/wells/StandardWellPrimaryVariables.hpp

@@ -102,7 +102,7 @@ public:
 
     //! \brief Copy values from well state.
     void update(const WellState& well_state,
-                const bool zero_rate_target,
+                const bool stop_or_zero_rate_target,
                 DeferredLogger& deferred_logger);
 
     //! \brief Copy polymer molecular weigt values from well state.
@@ -110,7 +110,7 @@ public:
 
     //! \brief Update values from newton update vector.
     void updateNewton(const BVectorWell& dwells,
-                      const bool zero_rate_target,
+                      const bool stop_or_zero_rate_target,
                       const double dFLimit,
                       const double dBHPLimit);
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -1000,10 +1000,10 @@ namespace Opm
     {
         if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
-        const bool zero_rate_target = this->wellUnderZeroProductionRateControl(summary_state, well_state);
+        const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
-        updatePrimaryVariablesNewton(dwells, zero_rate_target, deferred_logger);
+        updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
 
-        updateWellStateFromPrimaryVariables(well_state, deferred_logger);
+        updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, deferred_logger);
         Base::calculateReservoirRates(well_state.well(this->index_of_well_));
     }
 
@@ -1015,12 +1015,12 @@ namespace Opm
     void
     StandardWell<TypeTag>::
     updatePrimaryVariablesNewton(const BVectorWell& dwells,
-                                 const bool zero_rate_target,
+                                 const bool stop_or_zero_rate_target,
                                  DeferredLogger& deferred_logger)
     {
         const double dFLimit = this->param_.dwell_fraction_max_;
         const double dBHPLimit = this->param_.dbhp_max_rel_;
-        this->primary_variables_.updateNewton(dwells, zero_rate_target, dFLimit, dBHPLimit);
+        this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit);
 
         // for the water velocity and skin pressure
         if constexpr (Base::has_polymermw) {
@@ -1037,9 +1037,11 @@ namespace Opm
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
-    updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger) const
+    updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
+                                        WellState& well_state,
+                                        DeferredLogger& deferred_logger) const
     {
-        this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, deferred_logger);
+        this->StdWellEval::updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, deferred_logger);
 
         // other primary variables related to polymer injectivity study
         if constexpr (Base::has_polymermw) {
@@ -1982,8 +1984,8 @@ namespace Opm
     {
         if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
-        const bool zero_rate_target = this->wellUnderZeroProductionRateControl(summary_state, well_state);
+        const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
-        this->primary_variables_.update(well_state, zero_rate_target, deferred_logger);
+        this->primary_variables_.update(well_state, stop_or_zero_rate_target, deferred_logger);
 
         // other primary variables related to polymer injection
         if constexpr (Base::has_polymermw) {

opm/simulators/wells/WellBhpThpCalculator.cpp

@@ -224,6 +224,7 @@ computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>&
 }
 
 void WellBhpThpCalculator::updateThp(const double rho,
+                                     const bool stop_or_zero_rate_target,
                                      const std::function<double()>& alq_value,
                                      const std::array<unsigned,3>& active,
                                      WellState& well_state,
@@ -235,7 +236,7 @@ void WellBhpThpCalculator::updateThp(const double rho,
     auto& ws = well_state.well(well_.indexOfWell());
 
     // When there is no vaild VFP table provided, we set the thp to be zero.
-    if (!well_.isVFPActive(deferred_logger) || well_.wellIsStopped()) {
+    if (!well_.isVFPActive(deferred_logger) || stop_or_zero_rate_target) {
         ws.thp = 0;
         return;
     }

opm/simulators/wells/WellBhpThpCalculator.hpp

@@ -82,6 +82,7 @@ public:
 
     //! \brief Update THP.
     void updateThp(const double rho,
+                   const bool stop_or_zero_rate_target,
                    const std::function<double()>& alq_value,
                    const std::array<unsigned,3>& active,
                    WellState& well_state,

opm/simulators/wells/WellHelpers.cpp

@@ -25,6 +25,7 @@
 
 #include <opm/common/OpmLog/OpmLog.hpp>
 
+#include <opm/input/eclipse/Schedule/Well/WellInjectionControls.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellProductionControls.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellEnums.hpp>
 
@@ -171,8 +172,8 @@ DenseMatrix transposeDenseDynMatrix(const DenseMatrix& M)
     return tmp;
 }
 
-bool rateControlWithZeroTarget(const WellProducerCMode& mode,
+bool rateControlWithZeroProdTarget(const WellProductionControls& controls,
-                               const WellProductionControls& controls)
+                                   const WellProducerCMode mode)
 {
     switch (mode) {
         case WellProducerCMode::ORAT:
@@ -198,6 +199,20 @@ bool rateControlWithZeroTarget(const WellProducerCMode& mode,
 }
 
 
+bool rateControlWithZeroInjTarget(const WellInjectionControls& controls,
+                                  const WellInjectorCMode mode)
+{
+    switch (mode) {
+        case WellInjectorCMode::RATE:
+            return controls.surface_rate == 0.0;
+        case WellInjectorCMode::RESV:
+            return controls.reservoir_rate == 0.0;
+        default:
+            return false;
+    }
+}
+
+
 template class ParallelStandardWellB<double>;
 
 template<int Dim> using Vec = Dune::BlockVector<Dune::FieldVector<double,Dim>>;

opm/simulators/wells/WellHelpers.hpp

@@ -31,8 +31,10 @@
 namespace Opm {
 
 class ParallelWellInfo;
+struct WellProductionControls;
+struct WellInjectionControls;
 enum class WellProducerCMode;
-class WellProductionControls;
+enum class WellInjectorCMode;
 
 namespace wellhelpers {
 
@@ -86,8 +88,12 @@ template <class DenseMatrix>
 DenseMatrix transposeDenseDynMatrix(const DenseMatrix& M);
 
 /// Helper to check whether the well is under zero production rate control
-bool rateControlWithZeroTarget(const WellProducerCMode& mode,
+bool rateControlWithZeroProdTarget(const WellProductionControls& controls,
-                               const WellProductionControls& controls);
+                                   WellProducerCMode mode);
+
+/// Helper to check whether the well is under zero injection rate control
+bool rateControlWithZeroInjTarget(const WellInjectionControls& controls,
+                                  WellInjectorCMode mode);
 
 } // namespace wellhelpers
 } // namespace Opm

opm/simulators/wells/WellInterface.hpp

@@ -360,8 +360,6 @@ protected:
 
     Eval getPerfCellPressure(const FluidState& fs) const;
 
-    bool wellUnderZeroProductionRateControl(const SummaryState& summary_state,
-                                            const WellState& well_state) const;
 };
 
 }

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -450,6 +450,29 @@ bool WellInterfaceGeneric::isPressureControlled(const WellState& well_state) con
     }
 }
 
+
+
+bool WellInterfaceGeneric::wellUnderZeroRateTarget(const SummaryState& summary_state,
+                                                   const WellState& well_state) const
+{
+    if (this->isProducer()) { // producers
+        const auto prod_controls = this->well_ecl_.productionControls(summary_state);
+        const auto prod_mode = well_state.well(this->indexOfWell()).production_cmode;
+        return wellhelpers::rateControlWithZeroProdTarget(prod_controls, prod_mode);
+    } else { // injectors
+        const auto inj_controls = this->well_ecl_.injectionControls(summary_state);
+        const auto inj_mode = well_state.well(this->indexOfWell()).injection_cmode;
+        return wellhelpers::rateControlWithZeroInjTarget(inj_controls, inj_mode);
+    }
+}
+
+bool WellInterfaceGeneric::stopppedOrZeroRateTarget(const SummaryState& summary_state,
+                                                    const WellState& well_state) const
+{
+    return (this->wellIsStopped() || this->wellUnderZeroRateTarget(summary_state, well_state));
+
+}
+
 double WellInterfaceGeneric::wmicrobes_() const
 {
     auto injectorType = this->well_ecl_.injectorType();

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -195,6 +195,8 @@ public:
 
     bool isPressureControlled(const WellState& well_state) const;
 
+    bool stopppedOrZeroRateTarget(const SummaryState& summary_state,
+                                  const WellState& well_state) const;
 protected:
     bool getAllowCrossFlow() const;
 
@@ -209,6 +211,9 @@ protected:
     int polymerInjTable_() const;
     int polymerWaterTable_() const;
 
+    bool wellUnderZeroRateTarget(const SummaryState& summary_state,
+                                 const WellState& well_state) const;
+
     // definition of the struct OperabilityStatus
     struct OperabilityStatus {
         bool isOperableAndSolvable() const {

opm/simulators/wells/WellInterface_impl.hpp

@@ -184,7 +184,8 @@ namespace Opm
                       const GroupState& group_state,
                       DeferredLogger& deferred_logger) /* const */
     {
-        if (this->wellIsStopped()) {
+        const auto& summary_state = ebos_simulator.vanguard().summaryState();
+        if (this->stopppedOrZeroRateTarget(summary_state, well_state)) {
             return false;
         }
 
@@ -1131,18 +1132,4 @@ namespace Opm
         }
     }
 
-
-    template<typename TypeTag>
-    bool WellInterface<TypeTag>::wellUnderZeroProductionRateControl(const SummaryState& summary_state,
-                                                                    const WellState& well_state) const
-    {
-        if (this->wellIsStopped()) return true;
-
-        bool zero_rate_target = false;
-        if (this->well_ecl_.isProducer()) {
-            const auto prod_controls = this->well_ecl_.productionControls(summary_state);
-            zero_rate_target = wellhelpers::rateControlWithZeroTarget(well_state.well(this->index_of_well_).production_cmode, prod_controls);
-        }
-        return zero_rate_target;
-    }
 } // namespace Opm