Merge pull request #1648 from GitPaean/welloperability_checking

[WIP] Welloperability checking

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -396,6 +396,10 @@ namespace Opm {
         initial_step_ = false;
     }
 
+
+
+
+
     template<typename TypeTag>
     std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr >
     BlackoilWellModel<TypeTag>::
@@ -429,21 +433,33 @@ namespace Opm {
 
                 const Well* well_ecl = wells_ecl_[index_well];
 
-                // well is closed due to economical reasons
-                if (wellTestState_.hasWell(well_name, WellTestConfig::Reason::ECONOMIC)) { 
-		    if( well_ecl->getAutomaticShutIn() ) {
-                        // shut wells are not added to the well container 
-                        well_state_.bhp()[w] = 0;
+                // A new WCON keywords can re-open a well that was closed/shut due to Physical limit
+                if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL ) ) {
+                    // TODO: more checking here, to make sure this standard more specific and complete
+                    // maybe there is some WCON keywords will not open the well
+                    if (well_state_.effectiveEventsOccurred(w) ) {
+                        wellTestState_.openWell(well_name);
+                    }
+                }
+
+                // TODO: should we do this for all kinds of closing reasons?
+                // something like wellTestState_.hasWell(well_name)?
+                if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::ECONOMIC) ||
+                     wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL) ) {
+                    if( well_ecl->getAutomaticShutIn() ) {
+                        // shut wells are not added to the well container
+                        // TODO: make a function from well_state side to handle the following
+                        well_state_.thp()[w] = 0.;
+                        well_state_.bhp()[w] = 0.;
                         const int np = numPhases();
                         for (int p = 0; p < np; ++p) {
-                            well_state_.wellRates()[np * w + p] = 0;
-			}
-			continue;
-                    }
-		    else {
-		        // close wells are added to the container but marked as closed
-		        struct WellControls* well_controls = wells()->ctrls[w];
-		        well_controls_stop_well(well_controls);
+                            well_state_.wellRates()[np * w + p] = 0.;
+                        }
+                        continue;
+                    } else {
+                        // close wells are added to the container but marked as closed
+                        struct WellControls* well_controls = wells()->ctrls[w];
+                        well_controls_stop_well(well_controls);
                     }
                 }
 
@@ -803,7 +819,9 @@ namespace Opm {
         // Get global (from all processes) convergence report.
         ConvergenceReport local_report;
         for (const auto& well : well_container_) {
-            local_report += well->getWellConvergence(B_avg);
+            if (well->isOperable() ) {
+                local_report += well->getWellConvergence(B_avg);
+            }
         }
         ConvergenceReport report = gatherConvergenceReport(local_report);
 
@@ -828,9 +846,10 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     calculateExplicitQuantities() const
     {
-         for (auto& well : well_container_) {
-             well->calculateExplicitQuantities(ebosSimulator_, well_state_);
-         }
+        // TODO: checking isOperable() ?
+        for (auto& well : well_container_) {
+            well->calculateExplicitQuantities(ebosSimulator_, well_state_);
+        }
     }
 
 
@@ -934,6 +953,10 @@ namespace Opm {
         // process group control related
         prepareGroupControl();
 
+        for (const auto& well : well_container_) {
+            well->checkWellOperability(ebosSimulator_, well_state_);
+        }
+
         // since the controls are all updated, we should update well_state accordingly
         for (const auto& well : well_container_) {
             const int w = well->indexOfWell();
@@ -941,16 +964,22 @@ namespace Opm {
             const int control = well_controls_get_current(wc);
             well_state_.currentControls()[w] = control;
 
+            if (!well->isOperable() ) continue;
+
             if (well_state_.effectiveEventsOccurred(w) ) {
                 well->updateWellStateWithTarget(ebosSimulator_, well_state_);
             }
 
             // there is no new well control change input within a report step,
             // so next time step, the well does not consider to have effective events anymore
+            // TODO: if we can know whether this is the first time step within the report step,
+            // we do not need to set it to false
+            // TODO: we should do this at the end of the time step in case we will need it within
+            // this time step somewhere
             if (well_state_.effectiveEventsOccurred(w) ) {
                 well_state_.setEffectiveEventsOccurred(w, false);
             }
-        }  // end of for (int w = 0; w < nw; ++w)
+        }  // end of for (const auto& well : well_container_)
 
         updatePrimaryVariables();
     }

opm/autodiff/MultisegmentWell.hpp

@@ -330,6 +330,11 @@ namespace Opm
         // handling the overshooting and undershooting of the fractions
         void processFractions(const int seg) const;
 
+        // checking the operability of the well based on current reservoir condition
+        // it is not implemented for multisegment well yet
+        virtual void checkWellOperability(const Simulator& ebos_simulator,
+                                          const WellState& well_state) override;
+
         void updateWellStateFromPrimaryVariables(WellState& well_state) const;
 
         bool frictionalPressureLossConsidered() const;

opm/autodiff/MultisegmentWell_impl.hpp

@@ -1642,6 +1642,21 @@ namespace Opm
 
 
 
+    template <typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    checkWellOperability(const Simulator& /* ebos_simulator */,
+                         const WellState& /* well_state */)
+    {
+        const std::string msg = "Support of well operatability checking for mutlisegment wells is not implemented "
+                                "yet, checkWellOperability() for " + name() + " will do nothing";
+        OpmLog::warning("NO_OPERATABILITY_CHECKING_MS_WELLS", msg);
+    }
+
+
+
+
+
     template <typename TypeTag>
     void
     MultisegmentWell<TypeTag>::

opm/autodiff/StandardWell.hpp

@@ -264,8 +264,6 @@ namespace Opm
 
         EvalWell extendEval(const Eval& in) const;
 
-        bool crossFlowAllowed(const Simulator& ebosSimulator) const;
-
         // xw = inv(D)*(rw - C*x)
         void recoverSolutionWell(const BVector& x, BVectorWell& xw) const;
 
@@ -350,6 +348,17 @@ namespace Opm
         // updating the inflow based on the current reservoir condition
         void updateIPR(const Simulator& ebos_simulator) const;
 
+        // check whether the well is operable under the current reservoir condition
+        // mostly related to BHP limit and THP limit
+        virtual void checkWellOperability(const Simulator& ebos_simulator,
+                                          const WellState& well_state) override;
+
+        // check whether the well is operable under BHP limit with current reservoir condition
+        void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator);
+
+        // check whether the well is operable under THP limit with current reservoir condition
+        void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator);
+
         // update WellState based on IPR and associated VFP table
         void updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator,
                                              WellState& well_state) const;
@@ -357,6 +366,21 @@ namespace Opm
         void updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator,
                                                      WellState& well_state) const;
 
+        // for a well, when all drawdown are in the wrong direction, then this well will not
+        // be able to produce/inject .
+        bool allDrawDownWrongDirection(const Simulator& ebos_simulator) const;
+
+        // whether the well can produce / inject based on the current well state (bhp)
+        bool canProduceInjectWithCurrentBhp(const Simulator& ebos_simulator,
+                                            const WellState& well_state);
+
+        // turn on crossflow to avoid singular well equations
+        // when the well is banned from cross-flow and the BHP is not properly initialized,
+        // we turn on crossflow to avoid singular well equations. It can result in wrong-signed
+        // well rates, it can cause problem for THP calculation
+        // TODO: looking for better alternative to avoid wrong-signed well rates
+        bool openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const;
+
         // calculate the BHP from THP target based on IPR
         // TODO: we need to check the operablility here first, if not operable, then maybe there is
         // no point to do this

opm/autodiff/StandardWell_impl.hpp

@@ -440,9 +440,9 @@ namespace Opm
                    WellState& well_state)
     {
 
-        const Opm::SummaryConfig& summaryConfig = ebosSimulator.vanguard().summaryConfig();
+        checkWellOperability(ebosSimulator, well_state);
 
-        const int np = number_of_phases_;
+        if (!this->isOperable()) return;
 
         // clear all entries
         duneB_ = 0.0;
@@ -453,7 +453,7 @@ namespace Opm
         // TODO: it probably can be static member for StandardWell
         const double volume = 0.002831684659200; // 0.1 cu ft;
 
-        const bool allow_cf = crossFlowAllowed(ebosSimulator);
+        const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
 
         const EvalWell& bhp = getBhp();
 
@@ -461,6 +461,7 @@ namespace Opm
         well_state.wellVaporizedOilRates()[index_of_well_] = 0.;
         well_state.wellDissolvedGasRates()[index_of_well_] = 0.;
 
+        const int np = number_of_phases_;
         for (int p = 0; p < np; ++p) {
             well_state.productivityIndex()[np*index_of_well_ + p] = 0.;
         }
@@ -588,7 +589,8 @@ namespace Opm
             well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf];
 
             // Compute Productivity index if asked for
-             const auto& pu = phaseUsage();
+            const auto& pu = phaseUsage();
+            const Opm::SummaryConfig& summaryConfig = ebosSimulator.vanguard().summaryConfig();
             for (int p = 0; p < np; ++p) {
                 if ( (pu.phase_pos[Water] == p && (summaryConfig.hasSummaryKey("WPIW:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name())))
                         || (pu.phase_pos[Oil] == p && (summaryConfig.hasSummaryKey("WPIO:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name())))
@@ -747,45 +749,6 @@ namespace Opm
 
 
 
-    template<typename TypeTag>
-    bool
-    StandardWell<TypeTag>::
-    crossFlowAllowed(const Simulator& ebosSimulator) const
-    {
-        if (getAllowCrossFlow()) {
-            return true;
-        }
-
-        // TODO: investigate the justification of the following situation
-
-        // check for special case where all perforations have cross flow
-        // then the wells must allow for cross flow
-        for (int perf = 0; perf < number_of_perforations_; ++perf) {
-            const int cell_idx = well_cells_[perf];
-            const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
-            const auto& fs = intQuants.fluidState();
-            const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
-            const EvalWell& bhp = getBhp();
-
-            // Pressure drawdown (also used to determine direction of flow)
-            const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
-            const EvalWell drawdown = pressure - well_pressure;
-
-            if (drawdown.value() < 0 && well_type_ == INJECTOR)  {
-                return false;
-            }
-
-            if (drawdown.value() > 0 && well_type_ == PRODUCER)  {
-                return false;
-            }
-        }
-        return true;
-    }
-
-
-
-
-
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
@@ -860,6 +823,8 @@ namespace Opm
     updateWellState(const BVectorWell& dwells,
                     WellState& well_state) const
     {
+        if (!this->isOperable()) return;
+
         updatePrimaryVariablesNewton(dwells, well_state);
 
         updateWellStateFromPrimaryVariables(well_state);
@@ -1144,13 +1109,21 @@ namespace Opm
             // TODO: we should address this in a function updateWellStateWithBHPTarget.
             // TODO: however, the reason that this one minght not be that critical with
             // TODO: the effects remaining to be investigated.
-
             break;
+
         case THP: {
-            // TODO: adding the checking for the operability
-            // TODO: should we do updateIPR before this or within the related functions
-            updateIPR(ebos_simulator);
-            updateWellStateWithTHPTargetIPR(ebos_simulator, well_state);
+            assert(this->isOperable() );
+
+            // when a well can not work under THP target, it switches to BHP control
+            if (this->operability_status_.isOperableUnderTHPLimit() ) {
+                updateWellStateWithTHPTargetIPR(ebos_simulator, well_state);
+            } else { // go to BHP limit
+                assert(this->operability_status_.isOperableUnderBHPLimit() );
+
+                OpmLog::info("well " + name() + " can not work with THP target, switching to BHP control");
+
+                well_state.bhp()[well_index] = mostStrictBhpFromBhpLimits();
+            }
             break;
         }
 
@@ -1267,7 +1240,8 @@ namespace Opm
             // it should not be negative anyway. If it is negative, we might need to re-formulate
             // to taking into consideration the crossflow here.
             if (pressure_diff <= 0.) {
-                OpmLog::warning("NON_POSITIVE_DRAWDOWN_IPR", "non-positive drawdown found when updateIPR for well " + name());
+                OpmLog::warning("NON_POSITIVE_DRAWDOWN_IPR",
+                                "non-positive drawdown found when updateIPR for well " + name());
             }
 
             // the well index associated with the connection
@@ -1316,6 +1290,221 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    checkWellOperability(const Simulator& ebos_simulator,
+                         const WellState& well_state)
+    {
+        // focusing on PRODUCER for now
+        if (well_type_ == INJECTOR) {
+            return;
+        }
+
+        if (!this->underPredictionMode() ) {
+            return;
+        }
+
+        const bool old_well_operable = this->operability_status_.isOperable();
+
+        this->operability_status_.reset();
+
+        updateIPR(ebos_simulator);
+
+        // checking the BHP limit related
+        checkOperabilityUnderBHPLimitProducer(ebos_simulator);
+
+        // checking whether the well can operate under the THP constraints.
+        if (this->wellHasTHPConstraints()) {
+            this->operability_status_.has_thp_constaint = true;
+            checkOperabilityUnderTHPLimitProducer(ebos_simulator);
+            this->operability_status_.can_produce_inject_with_current_bhp =
+                canProduceInjectWithCurrentBhp(ebos_simulator, well_state);
+        }
+
+        const bool well_operable = this->operability_status_.isOperable();
+
+        if (!well_operable && old_well_operable) {
+            OpmLog::info(" well " + name() + " gets SHUT during iteration ");
+        } else if (well_operable && !old_well_operable) {
+            OpmLog::info(" well " + name() + " gets REVIVED during iteration ");
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator)
+    {
+        const double bhp_limit = mostStrictBhpFromBhpLimits();
+        // Crude but works: default is one atmosphere.
+        // TODO: a better way to detect whether the BHP is defaulted or not
+        const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
+        if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints() ) {
+            // if the BHP limit is not defaulted or the well does not have a THP limit
+            // we need to check the BHP limit
+
+            for (int p = 0; p < number_of_phases_; ++p) {
+                const double temp = ipr_a_[p] - ipr_b_[p] * bhp_limit;
+                if (temp < 0.) {
+                    this->operability_status_.operable_under_only_bhp_limit = false;
+                    break;
+                }
+            }
+
+            // checking whether running under BHP limit will violate THP limit
+            if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints()) {
+                // option 1: calculate well rates based on the BHP limit.
+                // option 2: stick with the above IPR curve
+                // we use IPR here
+                std::vector<double> well_rates_bhp_limit;
+                computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit);
+
+                const double thp = calculateThpFromBhp(well_rates_bhp_limit, bhp_limit);
+                const double thp_limit = this->getTHPConstraint();
+
+                if (thp < thp_limit) {
+                    this->operability_status_.obey_thp_limit_under_bhp_limit = false;
+                }
+            }
+        } else {
+            // defaulted BHP and there is a THP constraint
+            // default BHP limit is about 1 atm.
+            // when applied the hydrostatic pressure correction dp,
+            // most likely we get a negative value (bhp + dp)to search in the VFP table,
+            // which is not desirable.
+            // we assume we can operate under defaulted BHP limit and will violate the THP limit
+            // when operating under defaulted BHP limit.
+            this->operability_status_.operable_under_only_bhp_limit = true;
+            this->operability_status_.obey_thp_limit_under_bhp_limit = false;
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator)
+    {
+        const double obtain_bhp =  calculateBHPWithTHPTargetIPR();
+
+        if (obtain_bhp > 0.) {
+            this->operability_status_.can_obtain_bhp_with_thp_limit = true;
+
+            const double  bhp_limit = mostStrictBhpFromBhpLimits();
+            this->operability_status_.obey_bhp_limit_with_thp_limit = (obtain_bhp >= bhp_limit);
+
+            const double thp_limit = this->getTHPConstraint();
+            if (obtain_bhp < thp_limit) {
+                const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(obtain_bhp, unit::barsa))
+                                        + " bars is SMALLER than thp limit "
+                                        + std::to_string(unit::convert::to(thp_limit, unit::barsa))
+                                        + " bars as a producer for well " + name();
+                OpmLog::debug(msg);
+            }
+        } else {
+            this->operability_status_.can_obtain_bhp_with_thp_limit = false;
+            const double thp_limit = this->getTHPConstraint();
+            OpmLog::debug(" COULD NOT find bhp value under thp_limit "
+                          + std::to_string(unit::convert::to(thp_limit, unit::barsa))
+                          + " bars for well " + name() + ", the well might need to be closed ");
+            this->operability_status_.obey_bhp_limit_with_thp_limit = false;
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    bool
+    StandardWell<TypeTag>::
+    allDrawDownWrongDirection(const Simulator& ebos_simulator) const
+    {
+        bool all_drawdown_wrong_direction = true;
+
+        for (int perf = 0; perf < number_of_perforations_; ++perf) {
+            const int cell_idx = well_cells_[perf];
+            const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
+            const auto& fs = intQuants.fluidState();
+
+            const double pressure = (fs.pressure(FluidSystem::oilPhaseIdx)).value();
+            const double bhp = getBhp().value();
+
+            // Pressure drawdown (also used to determine direction of flow)
+            const double well_pressure = bhp + perf_pressure_diffs_[perf];
+            const double drawdown = pressure - well_pressure;
+
+            // for now, if there is one perforation can produce/inject in the correct
+            // direction, we consider this well can still produce/inject.
+            // TODO: it can be more complicated than this to cause wrong-signed rates
+            if ( (drawdown < 0. && well_type_ == INJECTOR) ||
+                 (drawdown > 0. && well_type_ == PRODUCER) )  {
+                all_drawdown_wrong_direction = false;
+                break;
+            }
+        }
+
+        return all_drawdown_wrong_direction;
+    }
+
+
+
+
+    template<typename TypeTag>
+    bool
+    StandardWell<TypeTag>::
+    canProduceInjectWithCurrentBhp(const Simulator& ebos_simulator,
+                                   const WellState& well_state)
+    {
+        const double bhp = well_state.bhp()[index_of_well_];
+        std::vector<double> well_rates;
+        computeWellRatesWithBhp(ebos_simulator, bhp, well_rates);
+
+        const double sign = (well_type_ == PRODUCER) ? -1. : 1.;
+        const double threshold = sign * std::numeric_limits<double>::min();
+
+        bool can_produce_inject = false;
+        for (const auto value : well_rates) {
+            if (well_type_ == PRODUCER && value < threshold) {
+                can_produce_inject = true;
+                break;
+            } else if (well_type_ == INJECTOR && value > threshold) {
+                can_produce_inject = true;
+                break;
+            }
+        }
+
+        if (!can_produce_inject) {
+            OpmLog::debug(" well " + name() + " CANNOT produce or inejct ");
+        }
+
+        return can_produce_inject;
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    bool
+    StandardWell<TypeTag>::
+    openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const
+    {
+        return !getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
@@ -1349,6 +1538,8 @@ namespace Opm
 
         const double bhp = calculateBHPWithTHPTargetIPR();
 
+        assert(bhp > 0.0);
+
         well_state.bhp()[index_of_well_] = bhp;
 
         // TODO: explicit quantities are always tricky for this type of situation
@@ -1401,9 +1592,6 @@ namespace Opm
         const double obtain_bhp = vfp_properties_->getProd()->calculateBhpWithTHPTarget(ipr_a_, ipr_b_,
                                              bhp_limit, thp_table_id, thp_target, alq, dp);
 
-        // we should have made sure that this well should be operable under THP limit now
-        assert(obtain_bhp > 0.);
-
         return obtain_bhp;
     }
 
@@ -1726,7 +1914,7 @@ namespace Opm
         switch(well_controls_get_current_type(well_controls_)) {
             case THP:
                 type = CR::WellFailure::Type::ControlTHP;
-                control_tolerance = 1.e3; // 0.01 bar
+                control_tolerance = 1.e4; // 0.1 bar
                 break;
             case BHP:  // pressure type of control
                 type = CR::WellFailure::Type::ControlBHP;
@@ -1816,6 +2004,8 @@ namespace Opm
     StandardWell<TypeTag>::
     solveEqAndUpdateWellState(WellState& well_state)
     {
+        if (!this->isOperable()) return;
+
         // We assemble the well equations, then we check the convergence,
         // which is why we do not put the assembleWellEq here.
         BVectorWell dx_well(1);
@@ -1861,6 +2051,8 @@ namespace Opm
     StandardWell<TypeTag>::
     apply(const BVector& x, BVector& Ax) const
     {
+        if (!this->isOperable()) return;
+
         if ( param_.matrix_add_well_contributions_ )
         {
             // Contributions are already in the matrix itself
@@ -1889,6 +2081,8 @@ namespace Opm
     StandardWell<TypeTag>::
     apply(BVector& r) const
     {
+        if (!this->isOperable()) return;
+
         assert( invDrw_.size() == invDuneD_.N() );
 
         // invDrw_ = invDuneD_ * resWell_
@@ -1906,6 +2100,8 @@ namespace Opm
     StandardWell<TypeTag>::
     recoverSolutionWell(const BVector& x, BVectorWell& xw) const
     {
+        if (!this->isOperable()) return;
+
         BVectorWell resWell = resWell_;
         // resWell = resWell - B * x
         duneB_.mmv(x, resWell);
@@ -1923,6 +2119,8 @@ namespace Opm
     recoverWellSolutionAndUpdateWellState(const BVector& x,
                                           WellState& well_state) const
     {
+        if (!this->isOperable()) return;
+
         BVectorWell xw(1);
         recoverSolutionWell(x, xw);
         updateWellState(xw, well_state);
@@ -1942,7 +2140,7 @@ namespace Opm
         const int np = number_of_phases_;
         well_flux.resize(np, 0.0);
 
-        const bool allow_cf = crossFlowAllowed(ebosSimulator);
+        const bool allow_cf = getAllowCrossFlow();
 
         for (int perf = 0; perf < number_of_perforations_; ++perf) {
             const int cell_idx = well_cells_[perf];
@@ -2127,6 +2325,8 @@ namespace Opm
     StandardWell<TypeTag>::
     updatePrimaryVariables(const WellState& well_state) const
     {
+        if (!this->isOperable()) return;
+
         const int well_index = index_of_well_;
         const int np = number_of_phases_;
 
@@ -2337,7 +2537,8 @@ namespace Opm
                 return;
             }
             // compute the well water velocity with out shear effects.
-            const bool allow_cf = crossFlowAllowed(ebos_simulator);
+            // TODO: do we need to turn on crossflow here?
+            const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
             const EvalWell& bhp = getBhp();
             std::vector<EvalWell> cq_s(num_components_,0.0);
             double perf_dis_gas_rate = 0.;

opm/autodiff/VFPProdProperties.cpp

@@ -229,10 +229,10 @@ calculateBhpWithTHPTarget(const std::vector<double>& ipr_a,
                                                                     detail::RateBhpPair{flo_bhp_limit, bhp_limit} };
 
     double obtain_bhp = 0.;
-    const bool obtain_solution_with_thp_limit = detail::findIntersectionForBhp(ratebhp_samples, ratebhp_twopoints_ipr, obtain_bhp);
+    const bool can_obtain_bhp_with_thp_limit = detail::findIntersectionForBhp(ratebhp_samples, ratebhp_twopoints_ipr, obtain_bhp);
 
-    // \Note: assuming not that negative BHP does not make sense
-    if (obtain_solution_with_thp_limit && obtain_bhp > 0.) {
+    // \Note: assuming that negative BHP does not make sense
+    if (can_obtain_bhp_with_thp_limit && obtain_bhp > 0.) {
         // getting too high bhp that might cause negative rates (rates in the undesired direction)
         if (obtain_bhp >= bhp_safe_limit) {
             const std::string msg (" We are getting a too high BHP value from the THP constraint, which may "
@@ -252,5 +252,4 @@ calculateBhpWithTHPTarget(const std::vector<double>& ipr_a,
 }
 
 
-
 }

opm/autodiff/WellInterface.hpp

@@ -223,6 +223,10 @@ namespace Opm
 
         void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
 
+        virtual void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state) = 0;
+
+        // whether the well is operable
+        bool isOperable() const;
 
     protected:
 
@@ -350,6 +354,10 @@ namespace Opm
         // whether a well is specified with a non-zero and valid VFP table number
         bool isVFPActive() const;
 
+        struct OperabilityStatus;
+
+        OperabilityStatus operability_status_;
+
         void wellTestingEconomic(Simulator& simulator, const std::vector<double>& B_avg,
                                  const double simulation_time, const int report_step, const bool terminal_output,
                                  const WellState& well_state, WellTestState& welltest_state);
@@ -359,6 +367,11 @@ namespace Opm
                                          const bool write_message_to_opmlog,
                                          WellTestState& well_test_state) const;
 
+        void updateWellTestStatePhysical(const WellState& well_state,
+                                         const double simulation_time,
+                                         const bool write_message_to_opmlog,
+                                         WellTestState& well_test_state) const;
+
         void solveWellForTesting(Simulator& ebosSimulator, WellState& well_state, const std::vector<double>& B_avg, bool terminal_output);
 
         void scaleProductivityIndex(const int perfIdx, double& productivity_index) const;
@@ -366,6 +379,62 @@ namespace Opm
 
     };
 
+
+
+
+    // definition of the struct OperabilityStatus
+    template<typename TypeTag>
+    struct
+    WellInterface<TypeTag>::
+    OperabilityStatus {
+        bool isOperable() const {
+            if (!operable_under_only_bhp_limit) {
+                return false;
+            } else {
+                return ( (isOperableUnderBHPLimit() || isOperableUnderTHPLimit()) &&
+                        !(has_thp_constaint && !can_produce_inject_with_current_bhp) );
+            }
+        }
+
+        bool isOperableUnderBHPLimit() const {
+            return operable_under_only_bhp_limit && obey_thp_limit_under_bhp_limit;
+        }
+
+        bool isOperableUnderTHPLimit() const {
+            return can_obtain_bhp_with_thp_limit && obey_bhp_limit_with_thp_limit;
+        }
+
+        void reset() {
+            operable_under_only_bhp_limit = true;
+            obey_thp_limit_under_bhp_limit = true;
+            can_obtain_bhp_with_thp_limit = true;
+            obey_bhp_limit_with_thp_limit = true;
+            can_produce_inject_with_current_bhp = true;
+            has_thp_constaint = false;
+        }
+
+        // whether the well can be operated under bhp limit
+        // without considering other limits.
+        // if it is false, then the well is not operable for sure.
+        bool operable_under_only_bhp_limit = true;
+        // if the well can be operated under bhp limit, will it obey(not violate)
+        // the thp limit when operated under bhp limit
+        bool obey_thp_limit_under_bhp_limit = true;
+        // whether the well operate under the thp limit only
+        bool can_obtain_bhp_with_thp_limit = true;
+        // whether the well obey bhp limit when operated under thp limit
+        bool obey_bhp_limit_with_thp_limit = true;
+
+        // TODO: the following criterion is based on the current state of
+        // the well, we consider it is a numerical criterion.
+        // at the moment, we only apply it with well has THP constraint.
+        // whether the well can produce / inject with the current bhp of the well
+        // it might be updated with other criterion with investigation with more cases.
+        bool can_produce_inject_with_current_bhp = true;
+        // whether the well has a THP constraint
+        bool has_thp_constaint = false;
+    };
+
 }
 
 #include "WellInterface_impl.hpp"

opm/autodiff/WellInterface_impl.hpp

@@ -457,8 +457,25 @@ namespace Opm
             if (wellhelpers::constraintBroken(
                     well_state.bhp(), well_state.thp(), well_state.wellRates(),
                     w, np, well_type_, wc, ctrl_index)) {
-                // ctrl_index will be the index of the broken constraint after the loop.
-                break;
+
+                // if the well can not work under THP / BHP control, we should not switch to THP / BHP control
+                const bool cannot_switch_to_bhp = well_controls_iget_type(wc, ctrl_index) == BHP && !operability_status_.isOperableUnderBHPLimit();
+                const bool cannot_switch_to_thp = well_controls_iget_type(wc, ctrl_index) == THP && !operability_status_.isOperableUnderTHPLimit();
+                const bool cannot_switch = cannot_switch_to_bhp || cannot_switch_to_thp;
+                if ( !cannot_switch ) {
+
+                    // ctrl_index will be the index of the broken constraint after the loop.
+                    break;
+                } else {
+                    // before we figure out to handle it, we give some debug information here
+                    if ( well_controls_iget_type(wc, ctrl_index) == BHP && !operability_status_.isOperableUnderBHPLimit() ) {
+                        OpmLog::debug("well " + name() + " breaks the BHP limit, while it is not operable under BHP limit");
+                    }
+
+                    if ( well_controls_iget_type(wc, ctrl_index) == THP && !operability_status_.isOperableUnderTHPLimit() ) {
+                        OpmLog::debug("well " + name() + " breaks the THP limit, while it is not operable under THP limit");
+                    }
+                }
             }
         }
 
@@ -719,6 +736,15 @@ namespace Opm
             return;
         }
 
+        // Based on current understanding, only under prediction mode, we need to shut well due to various
+        // reasons or limits. With more knowlage or testing cases later, this might need to be corrected.
+        if (!underPredictionMode() ) {
+            return;
+        }
+
+        // updating well test state based on physical (THP/BHP) limits.
+        updateWellTestStatePhysical(well_state, simulationTime, writeMessageToOPMLog, wellTestState);
+
         // updating well test state based on Economic limits.
         updateWellTestStateEconomic(well_state, simulationTime, writeMessageToOPMLog, wellTestState);
 
@@ -729,6 +755,30 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    updateWellTestStatePhysical(const WellState& well_state,
+                                const double simulation_time,
+                                const bool write_message_to_opmlog,
+                                WellTestState& well_test_state) const
+    {
+        if (!isOperable()) {
+            well_test_state.addClosedWell(name(), WellTestConfig::Reason::PHYSICAL, simulation_time);
+            if (write_message_to_opmlog) {
+                // TODO: considering auto shut in?
+                const std::string msg = "well " + name()
+                             + std::string(" will be shut as it can not operate under current reservoir condition");
+                OpmLog::info(msg);
+            }
+        }
+
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::
@@ -1207,4 +1257,15 @@ namespace Opm
     }
 
 
+
+
+
+    template<typename TypeTag>
+    bool
+    WellInterface<TypeTag>::
+    isOperable() const {
+        return operability_status_.isOperable();
+    }
+
+
 }