Add check for operability of MSW

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -343,8 +343,11 @@ namespace Opm {
             // do the initialization for all the wells
             // TODO: to see whether we can postpone of the intialization of the well containers to
             // optimize the usage of the following several member variables
+            std::vector< Scalar > B_avg(numComponents(), Scalar() );
+            computeAverageFormationFactor(B_avg);
+
             for (auto& well : well_container_) {
-                well->init(&phase_usage_, depth_, gravity_, local_num_cells_);
+                well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg);
             }
 
             // update the updated cell flag
@@ -435,7 +438,7 @@ namespace Opm {
                 WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx, deferred_logger);
 
                 // some preparation before the well can be used
-                well->init(&phase_usage_, depth_, gravity_, local_num_cells_);
+                well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg);
                 const Well& wellEcl = schedule().getWell(well_name, timeStepIdx);
                 double well_efficiency_factor = wellEcl.getEfficiencyFactor();
                 WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), timeStepIdx), schedule(), timeStepIdx, well_efficiency_factor);

opm/simulators/wells/MultisegmentWell.hpp

@@ -113,7 +113,8 @@ namespace Opm
         virtual void init(const PhaseUsage* phase_usage_arg,
                           const std::vector<double>& depth_arg,
                           const double gravity_arg,
-                          const int num_cells) override;
+                          const int num_cells,
+                          const std::vector< Scalar >& B_avg) override;
 
 
         virtual void initPrimaryVariablesEvaluation() const override;
@@ -230,6 +231,9 @@ namespace Opm
         using Base::perf_depth_;
         using Base::num_components_;
         using Base::connectionRates_;
+        using Base::ipr_a_;
+        using Base::ipr_b_;
+        using Base::changed_to_stopped_this_step_;
 
         // protected functions from the Base class
         using Base::phaseUsage;
@@ -239,7 +243,8 @@ namespace Opm
         using Base::getAllowCrossFlow;
         using Base::scalingFactor;
         using Base::wellIsStopped_;
-
+        using Base::updateWellOperability;
+        using Base::checkWellOperability;
 
         // TODO: trying to use the information from the Well opm-parser as much
         // as possible, it will possibly be re-implemented later for efficiency reason.
@@ -431,12 +436,6 @@ 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,
-                                          Opm::DeferredLogger& deferred_logger) override;
-
         void updateWellStateFromPrimaryVariables(WellState& well_state, Opm::DeferredLogger& deferred_logger) const;
 
         bool frictionalPressureLossConsidered() const;
@@ -458,10 +457,6 @@ namespace Opm
                                     WellState& well_state,
                                     Opm::DeferredLogger& deferred_logger) override;
 
-        virtual void wellTestingPhysical(const Simulator& simulator, const std::vector<double>& B_avg,
-                                         const double simulation_time, const int report_step,
-                                         WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) override;
-
         virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;
 
         EvalWell getSegmentSurfaceVolume(const Simulator& ebos_simulator, const int seg_idx) const;
@@ -517,6 +512,16 @@ namespace Opm
         void assembleValvePressureEq(const int seg, WellState& well_state) const;
 
         EvalWell pressureDropValve(const int seg) const;
+
+        // check whether the well is operable under BHP limit with current reservoir condition
+        virtual void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) override;
+
+        // check whether the well is operable under THP limit with current reservoir condition
+        virtual void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger) override;
+
+        // updating the inflow based on the current reservoir condition
+        virtual void updateIPR(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) const override;
+
     };
 
 }

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -135,9 +135,10 @@ namespace Opm
     init(const PhaseUsage* phase_usage_arg,
          const std::vector<double>& depth_arg,
          const double gravity_arg,
-         const int num_cells)
+         const int num_cells,
+         const std::vector< Scalar >& B_avg)
     {
-        Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells);
+        Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg);
 
         // TODO: for StandardWell, we need to update the perf depth here using depth_arg.
         // for MultisegmentWell, it is much more complicated.
@@ -264,6 +265,8 @@ namespace Opm
                    Opm::DeferredLogger& deferred_logger)
     {
 
+        checkWellOperability(ebosSimulator, well_state, deferred_logger);
+
         const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
         if (use_inner_iterations) {
             this->iterateWellEquations(ebosSimulator, B_avg, dt, well_state, deferred_logger);
@@ -624,6 +627,8 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     apply(const BVector& x, BVector& Ax) const
     {
+        if (!this->isOperable() && !this->wellIsStopped()) return;
+
         if ( param_.matrix_add_well_contributions_ )
         {
             // Contributions are already in the matrix itself
@@ -649,6 +654,8 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     apply(BVector& r) const
     {
+        if (!this->isOperable() && !this->wellIsStopped()) return;
+
         // invDrw_ = duneD^-1 * resWell_
         const BVectorWell invDrw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_);
         // r = r - duneC_^T * invDrw
@@ -726,6 +733,8 @@ namespace Opm
                                           WellState& well_state,
                                           Opm::DeferredLogger& deferred_logger) const
     {
+        if (!this->isOperable() && !this->wellIsStopped()) return;
+
         BVectorWell xw(1);
         recoverSolutionWell(x, xw);
         updateWellState(xw, well_state, deferred_logger);
@@ -820,7 +829,6 @@ namespace Opm
 
 
 
-
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
@@ -856,6 +864,7 @@ namespace Opm
             well_state_copy.currentProductionControls()[index_of_well_] = Well::ProducerCMode::BHP;
         }
         well_state_copy.bhp()[well_copy.index_of_well_] = bhp;
+
         well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
         const double dt = ebosSimulator.timeStepSize();
         // iterate to get a solution at the given bhp.
@@ -932,6 +941,7 @@ namespace Opm
     {
         // TODO: to test using rate conversion coefficients to see if it will be better than
         // this default one
+        if (!this->isOperable() && !this->wellIsStopped()) return;
 
         const Well& well = Base::wellEcl();
 
@@ -1005,6 +1015,8 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     recoverSolutionWell(const BVector& x, BVectorWell& xw) const
     {
+        if (!this->isOperable() && !this->wellIsStopped()) return;
+
         BVectorWell resWell = resWell_;
         // resWell = resWell - B * x
         duneB_.mmv(x, resWell);
@@ -1021,6 +1033,8 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     solveEqAndUpdateWellState(WellState& well_state, Opm::DeferredLogger& deferred_logger)
     {
+        if (!this->isOperable() && !this->wellIsStopped()) return;
+
         // We assemble the well equations, then we check the convergence,
         // which is why we do not put the assembleWellEq here.
         const BVectorWell dx_well = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_);
@@ -1113,6 +1127,8 @@ namespace Opm
                     Opm::DeferredLogger& deferred_logger,
                     const double relaxation_factor) const
     {
+        if (!this->isOperable() && !this->wellIsStopped()) return;
+
         const double dFLimit = param_.dwell_fraction_max_;
         const double max_pressure_change = param_.max_pressure_change_ms_wells_;
         const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;
@@ -2374,33 +2390,199 @@ namespace Opm
     }
 
 
-
+    template<typename TypeTag>
-
-
-    template <typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    checkWellOperability(const Simulator& /* ebos_simulator */,
+    checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger)
-                         const WellState& /* well_state */,
-                         Opm::DeferredLogger& deferred_logger)
     {
-        const bool checkOperability = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheck);
+        const auto& summaryState = ebos_simulator.vanguard().summaryState();
-        if (!checkOperability) {
+        const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
-            return;
+        // 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(summaryState) ) {
+            // if the BHP limit is not defaulted or the well does not have a THP limit
+            // we need to check the BHP limit
 
-        // focusing on PRODUCER for now
+            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(summaryState)) {
+                // 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, Base::B_avg_, bhp_limit, well_rates_bhp_limit, deferred_logger);
+
+                const double thp = calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, deferred_logger);
+
+                const double thp_limit = this->getTHPConstraint(summaryState);
+
+                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
+    MultisegmentWell<TypeTag>::
+    updateIPR(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) const
+    {
+        // TODO: not handling solvent related here for now
+
+        // TODO: it only handles the producers for now
+        // the formular for the injectors are not formulated yet
         if (this->isInjector()) {
             return;
         }
 
-        if (!this->underPredictionMode() ) {
+        // initialize all the values to be zero to begin with
-            return;
+        std::fill(ipr_a_.begin(), ipr_a_.end(), 0.);
-        }
+        std::fill(ipr_b_.begin(), ipr_b_.end(), 0.);
 
-        const std::string msg = "Support of well operability checking for multisegment wells is not implemented "
+        const int nseg = numberOfSegments();
-                                "yet, checkWellOperability() for " + name() + " will do nothing";
+        double seg_bhp_press_diff = 0;
-        deferred_logger.warning("NO_OPERATABILITY_CHECKING_MS_WELLS", msg);
+        double ref_depth = ref_depth_;
+        for (int seg = 0; seg < nseg; ++seg) {
+            // calculating the perforation rate for each perforation that belongs to this segment
+            const double segment_depth = segmentSet()[seg].depth();
+            const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth, segment_depth, segment_densities_[seg].value(), gravity_);
+            ref_depth = segment_depth;
+            seg_bhp_press_diff += dp;
+            for (const int perf : segment_perforations_[seg]) {
+            //std::vector<EvalWell> mob(num_components_, {numWellEq_ + numEq, 0.0});
+            std::vector<EvalWell> mob(num_components_, 0.0);
+
+            // TODO: mabye we should store the mobility somewhere, so that we only need to calculate it one per iteration
+            getMobility(ebos_simulator, perf, mob);
+
+            const int cell_idx = well_cells_[perf];
+            const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
+            const auto& fs = int_quantities.fluidState();
+            // the pressure of the reservoir grid block the well connection is in
+                    // pressure difference between the segment and the perforation
+            const double perf_seg_press_diff = gravity_ * segment_densities_[seg].value() * perforation_segment_depth_diffs_[perf];
+            // pressure difference between the perforation and the grid cell
+            const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf];
+            const double pressure_cell = fs.pressure(FluidSystem::oilPhaseIdx).value();
+
+            // calculating the b for the connection
+            std::vector<double> b_perf(num_components_);
+            for (size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
+                if (!FluidSystem::phaseIsActive(phase)) {
+                    continue;
+                }
+                const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
+                b_perf[comp_idx] = fs.invB(phase).value();
+            }
+
+            // the pressure difference between the connection and BHP
+            const double h_perf = cell_perf_press_diff + perf_seg_press_diff + seg_bhp_press_diff;
+            const double pressure_diff = pressure_cell - h_perf;
+
+            // Let us add a check, since the pressure is calculated based on zero value BHP
+            // 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.) {
+                deferred_logger.warning("NON_POSITIVE_DRAWDOWN_IPR",
+                                "non-positive drawdown found when updateIPR for well " + name());
+            }
+
+            // the well index associated with the connection
+            const double tw_perf = well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
+
+            // TODO: there might be some indices related problems here
+            // phases vs components
+            // ipr values for the perforation
+            std::vector<double> ipr_a_perf(ipr_a_.size());
+            std::vector<double> ipr_b_perf(ipr_b_.size());
+            for (int p = 0; p < number_of_phases_; ++p) {
+                const double tw_mob = tw_perf * mob[p].value() * b_perf[p];
+                ipr_a_perf[p] += tw_mob * pressure_diff;
+                ipr_b_perf[p] += tw_mob;
+            }
+
+            // we need to handle the rs and rv when both oil and gas are present
+            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
+                const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+                const double rs = (fs.Rs()).value();
+                const double rv = (fs.Rv()).value();
+
+                const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
+                const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
+
+                ipr_a_perf[gas_comp_idx] += dis_gas_a;
+                ipr_a_perf[oil_comp_idx] += vap_oil_a;
+
+                const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
+                const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
+
+                ipr_b_perf[gas_comp_idx] += dis_gas_b;
+                ipr_b_perf[oil_comp_idx] += vap_oil_b;
+            }
+
+            for (int p = 0; p < number_of_phases_; ++p) {
+                // TODO: double check the indices here
+                ipr_a_[ebosCompIdxToFlowCompIdx(p)] += ipr_a_perf[p];
+                ipr_b_[ebosCompIdxToFlowCompIdx(p)] += ipr_b_perf[p];
+            }
+            }
+        }
+    }
+
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger)
+    {
+        const auto& summaryState = ebos_simulator.vanguard().summaryState();
+        const auto obtain_bhp = computeBhpAtThpLimitProd(ebos_simulator, Base::B_avg_, summaryState, deferred_logger);
+
+        if (obtain_bhp) {
+            this->operability_status_.can_obtain_bhp_with_thp_limit = true;
+
+            const double  bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
+            this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
+
+            const double thp_limit = this->getTHPConstraint(summaryState);
+            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();
+                deferred_logger.debug(msg);
+            }
+        } else {
+            this->operability_status_.can_obtain_bhp_with_thp_limit = false;
+            this->operability_status_.obey_bhp_limit_with_thp_limit = false;
+            if (!this->wellIsStopped()) {
+                const double thp_limit = this->getTHPConstraint(summaryState);
+                deferred_logger.debug(" could not find bhp value at thp limit "
+                                      + std::to_string(unit::convert::to(thp_limit, unit::barsa))
+                                      + " bar for well " + name() + ", the well might need to be closed ");
+            }
+        }
     }
 
 
@@ -2503,6 +2685,8 @@ namespace Opm
                              WellState& well_state,
                              Opm::DeferredLogger& deferred_logger)
     {
+        if (!this->isOperable() && !this->wellIsStopped()) return true;
+
         const int max_iter_number = param_.max_inner_iter_ms_wells_;
         const WellState well_state0 = well_state;
         const std::vector<Scalar> residuals0 = getWellResiduals(B_avg);
@@ -2620,6 +2804,8 @@ namespace Opm
                                    Opm::DeferredLogger& deferred_logger)
     {
 
+        if (!this->isOperable() && !this->wellIsStopped()) return;
+
         // update the upwinding segments
         updateUpwindingSegments();
 
@@ -2847,20 +3033,6 @@ namespace Opm
     }
 
 
-    template<typename TypeTag>
-    void
-    MultisegmentWell<TypeTag>::
-    wellTestingPhysical(const Simulator& /* simulator */, const std::vector<double>& /* B_avg */,
-                        const double /* simulation_time */, const int /* report_step */,
-                        WellState& /* well_state */, WellTestState& /* welltest_state */, Opm::DeferredLogger& deferred_logger)
-    {
-        const std::string msg = "Support of well testing for physical limits for multisegment wells is not "
-                                "implemented yet, wellTestingPhysical() for " + name() + " will do nothing";
-        deferred_logger.warning("NO_WELLTESTPHYSICAL_CHECKING_MS_WELLS", msg);
-    }
-
-
-
 
 
     template<typename TypeTag>

opm/simulators/wells/StandardWell.hpp

@@ -174,7 +174,8 @@ namespace Opm
         virtual void init(const PhaseUsage* phase_usage_arg,
                           const std::vector<double>& depth_arg,
                           const double gravity_arg,
-                          const int num_cells) override;
+                          const int num_cells,
+                          const std::vector< Scalar >& B_avg) override;
 
 
         virtual void initPrimaryVariablesEvaluation() const override;
@@ -337,6 +338,8 @@ namespace Opm
         using Base::wfoam;
         using Base::scalingFactor;
         using Base::mostStrictBhpFromBhpLimits;
+        using Base::updateWellOperability;
+        using Base::checkWellOperability;
 
         // protected member variables from the Base class
         using Base::current_step_;
@@ -359,6 +362,9 @@ namespace Opm
         using Base::perf_rep_radius_;
         using Base::perf_length_;
         using Base::bore_diameters_;
+        using Base::ipr_a_;
+        using Base::ipr_b_;
+        using Base::changed_to_stopped_this_step_;
 
         using Base::wellIsStopped_;
 
@@ -397,14 +403,6 @@ namespace Opm
         // the saturations in the well bore under surface conditions at the beginning of the time step
         std::vector<double> F0_;
 
-        // the vectors used to describe the inflow performance relationship (IPR)
-        // Q = IPR_A - BHP * IPR_B
-        // TODO: it minght need to go to WellInterface, let us implement it in StandardWell first
-        // it is only updated and used for producers for now
-        mutable std::vector<double> ipr_a_;
-        mutable std::vector<double> ipr_b_;
-
-        bool changed_to_stopped_this_step_ = false;
         // Enable GLIFT debug mode. This will enable output of logging messages.
         bool glift_debug = false;
 
@@ -523,14 +521,6 @@ namespace Opm
         // handle the non reasonable fractions due to numerical overshoot
         void processFractions() const;
 
-        // updating the inflow based on the current reservoir condition
-        void updateIPR(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) const;
-
-        // update the operability status of 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,
-                                          Opm::DeferredLogger& deferred_logger) override;
 
         virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
                                                     const double dt,
@@ -553,17 +543,14 @@ namespace Opm
                                  EvalWell& cq_s_zfrac_effective,
                                  Opm::DeferredLogger& deferred_logger) const;
 
-        // check whether the well is operable under the current reservoir condition
-        // mostly related to BHP limit and THP limit
-        void updateWellOperability(const Simulator& ebos_simulator,
-                                   const WellState& well_state,
-                                   Opm::DeferredLogger& deferred_logger);
-
         // check whether the well is operable under BHP limit with current reservoir condition
-        void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
+        virtual void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) override;
 
         // check whether the well is operable under THP limit with current reservoir condition
-        void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger);
+        virtual void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger) override;
+
+        // updating the inflow based on the current reservoir condition
+        virtual void updateIPR(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) const override;
 
         // for a well, when all drawdown are in the wrong direction, then this well will not
         // be able to produce/inject .
@@ -594,11 +581,6 @@ namespace Opm
         static double relaxationFactorRate(const std::vector<double>& primary_variables,
                                            const BVectorWell& dwells);
 
-        virtual void wellTestingPhysical(const Simulator& simulator, const std::vector<double>& B_avg,
-                                         const double simulation_time, const int report_step,
-                                         WellState& well_state, WellTestState& welltest_state,
-                                         Opm::DeferredLogger& deferred_logger) override;
-
         // calculate the skin pressure based on water velocity, throughput and polymer concentration.
         // throughput is used to describe the formation damage during water/polymer injection.
         // calculated skin pressure will be applied to the drawdown during perforation rate calculation

opm/simulators/wells/StandardWell_impl.hpp

@@ -49,8 +49,6 @@ namespace Opm
     , perf_pressure_diffs_(number_of_perforations_)
     , parallelB_(duneB_, pw_info)
     , F0_(numWellConservationEq)
-    , ipr_a_(number_of_phases_)
-    , ipr_b_(number_of_phases_)
     {
         assert(num_components_ == numWellConservationEq);
 
@@ -69,9 +67,10 @@ namespace Opm
     init(const PhaseUsage* phase_usage_arg,
          const std::vector<double>& depth_arg,
          const double gravity_arg,
-         const int num_cells)
+         const int num_cells,
+         const std::vector< Scalar >& B_avg)
     {
-        Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells);
+        Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg);
 
         perf_depth_.resize(number_of_perforations_, 0.);
         for (int perf = 0; perf < number_of_perforations_; ++perf) {
@@ -1680,88 +1679,6 @@ namespace Opm
     }
 
 
-
-
-
-    template<typename TypeTag>
-    void
-    StandardWell<TypeTag>::
-    checkWellOperability(const Simulator& ebos_simulator,
-                         const WellState& well_state,
-                         Opm::DeferredLogger& deferred_logger)
-    {
-
-        const bool checkOperability = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheck);
-        if (!checkOperability) {
-            return;
-        }
-
-        // focusing on PRODUCER for now
-        if (this->isInjector()) {
-            return;
-        }
-
-        if (!this->underPredictionMode() ) {
-            return;
-        }
-
-        if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
-            return;
-        }
-
-        const bool old_well_operable = this->operability_status_.isOperable();
-
-        updateWellOperability(ebos_simulator, well_state, deferred_logger);
-
-        const bool well_operable = this->operability_status_.isOperable();
-
-        if (!well_operable && old_well_operable) {
-            if (well_ecl_.getAutomaticShutIn()) {
-                deferred_logger.info(" well " + name() + " gets SHUT during iteration ");
-            } else {
-                if (!this->wellIsStopped()) {
-                    deferred_logger.info(" well " + name() + " gets STOPPED during iteration ");
-                    this->stopWell();
-                    changed_to_stopped_this_step_ = true;
-                }
-            }
-        } else if (well_operable && !old_well_operable) {
-            deferred_logger.info(" well " + name() + " gets REVIVED during iteration ");
-            this->openWell();
-            changed_to_stopped_this_step_ = false;
-        }
-    }
-
-
-
-
-
-    template<typename TypeTag>
-    void
-    StandardWell<TypeTag>::
-    updateWellOperability(const Simulator& ebos_simulator,
-                          const WellState& well_state,
-                          Opm::DeferredLogger& deferred_logger)
-    {
-        this->operability_status_.reset();
-
-        updateIPR(ebos_simulator, deferred_logger);
-
-        // checking the BHP limit related
-        checkOperabilityUnderBHPLimitProducer(well_state, ebos_simulator, deferred_logger);
-
-        const auto& summaryState = ebos_simulator.vanguard().summaryState();
-
-        // checking whether the well can operate under the THP constraints.
-        if (this->wellHasTHPConstraints(summaryState)) {
-            checkOperabilityUnderTHPLimitProducer(ebos_simulator, well_state, deferred_logger);
-        }
-    }
-
-
-
-
-
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
@@ -3346,70 +3263,6 @@ namespace Opm
 
 
 
-
-    template<typename TypeTag>
-    void
-    StandardWell<TypeTag>::
-    wellTestingPhysical(const Simulator& ebos_simulator, const std::vector<double>& B_avg,
-                        const double /* simulation_time */, const int /* report_step */,
-                        WellState& well_state, WellTestState& welltest_state,
-                        Opm::DeferredLogger& deferred_logger)
-    {
-        deferred_logger.info(" well " + name() + " is being tested for physical limits");
-
-        // some most difficult things are the explicit quantities, since there is no information
-        // in the WellState to do a decent initialization
-
-        // TODO: Let us assume that the simulator is updated
-
-        // Let us try to do a normal simualtion running, to keep checking the operability status
-        // If the well is not operable during any of the time. It means it does not pass the physical
-        // limit test.
-
-        // create a copy of the well_state to use. If the operability checking is sucessful, we use this one
-        // to replace the original one
-        WellState well_state_copy = well_state;
-
-        // TODO: well state for this well is kind of all zero status
-        // we should be able to provide a better initialization
-        calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger);
-
-        updateWellOperability(ebos_simulator, well_state_copy, deferred_logger);
-
-        if ( !this->isOperable() ) {
-            const std::string msg = " well " + name() + " is not operable during well testing for physical reason";
-            deferred_logger.debug(msg);
-            return;
-        }
-
-        updateWellStateWithTarget(ebos_simulator, well_state_copy, deferred_logger);
-
-        calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger);
-
-        const double dt = ebos_simulator.timeStepSize();
-        const bool converged = this->iterateWellEquations(ebos_simulator, B_avg, dt, well_state_copy, deferred_logger);
-
-        if (!converged) {
-            const std::string msg = " well " + name() + " did not get converged during well testing for physical reason";
-            deferred_logger.debug(msg);
-            return;
-        }
-
-        if (this->isOperable() ) {
-            welltest_state.openWell(name(), WellTestConfig::PHYSICAL );
-            const std::string msg = " well " + name() + " is re-opened through well testing for physical reason";
-            deferred_logger.info(msg);
-            well_state = well_state_copy;
-        } else {
-            const std::string msg = " well " + name() + " is not operable during well testing for physical reason";
-            deferred_logger.debug(msg);
-        }
-    }
-
-
-
-
-
     template<typename TypeTag>
     typename StandardWell<TypeTag>::EvalWell
     StandardWell<TypeTag>::

opm/simulators/wells/WellInterface.hpp

@@ -158,7 +158,8 @@ namespace Opm
         virtual void init(const PhaseUsage* phase_usage_arg,
                           const std::vector<double>& depth_arg,
                           const double gravity_arg,
-                          const int num_cells);
+                          const int num_cells,
+                          const std::vector< Scalar >& B_avg);
 
         virtual void initPrimaryVariablesEvaluation() const = 0;
 
@@ -271,7 +272,14 @@ namespace Opm
 
         void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
 
-        virtual void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger) = 0;
+        void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger);
+
+        // check whether the well is operable under the current reservoir condition
+        // mostly related to BHP limit and THP limit
+        void updateWellOperability(const Simulator& ebos_simulator,
+                                   const WellState& well_state,
+                                   Opm::DeferredLogger& deferred_logger);
+
 
         // whether the well is operable
         bool isOperable() const;
@@ -411,6 +419,17 @@ namespace Opm
 
         std::optional<double> dynamic_thp_limit_;
 
+        std::vector< Scalar > B_avg_;
+
+        // the vectors used to describe the inflow performance relationship (IPR)
+        // Q = IPR_A - BHP * IPR_B
+        // TODO: it minght need to go to WellInterface, let us implement it in StandardWell first
+        // it is only updated and used for producers for now
+        mutable std::vector<double> ipr_a_;
+        mutable std::vector<double> ipr_b_;
+
+        bool changed_to_stopped_this_step_ = false;
+
         const PhaseUsage& phaseUsage() const;
 
         int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
@@ -476,13 +495,22 @@ namespace Opm
 
         OperabilityStatus operability_status_;
 
+        // check whether the well is operable under BHP limit with current reservoir condition
+        virtual void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) =0;
+
+        // check whether the well is operable under THP limit with current reservoir condition
+        virtual void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger) =0;
+
+        virtual void updateIPR(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) const=0;
+
+
         void wellTestingEconomic(const Simulator& simulator, const std::vector<double>& B_avg,
                                  const double simulation_time, const WellState& well_state,
                                  WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger);
 
-        virtual void wellTestingPhysical(const Simulator& simulator, const std::vector<double>& B_avg,
+        void wellTestingPhysical(const Simulator& simulator, const std::vector<double>& B_avg,
                                  const double simulation_time, const int report_step,
-                                         WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) = 0;
+                                 WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger);
 
 
         virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,

opm/simulators/wells/WellInterface_impl.hpp

@@ -51,6 +51,8 @@ namespace Opm
       , index_of_well_(index_of_well)
       , first_perf_(first_perf_index)
       , perf_data_(&perf_data)
+      , ipr_a_(number_of_phases_)
+      , ipr_b_(number_of_phases_)
     {
         assert(well.name()==pw_info.name());
         assert(std::is_sorted(perf_data.begin(), perf_data.end(),
@@ -122,10 +124,12 @@ namespace Opm
     init(const PhaseUsage* phase_usage_arg,
          const std::vector<double>& /* depth_arg */,
          const double gravity_arg,
-         const int /* num_cells */)
+         const int /* num_cells */,
+         const std::vector< Scalar >& B_avg)
     {
         phase_usage_ = phase_usage_arg;
         gravity_ = gravity_arg;
+        B_avg_ = B_avg;
     }
 
 
@@ -1409,6 +1413,141 @@ namespace Opm
     }
 
 
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    wellTestingPhysical(const Simulator& ebos_simulator, const std::vector<double>& B_avg,
+                        const double /* simulation_time */, const int /* report_step */,
+                        WellState& well_state, WellTestState& welltest_state,
+                        Opm::DeferredLogger& deferred_logger)
+    {
+        deferred_logger.info(" well " + name() + " is being tested for physical limits");
+
+        // some most difficult things are the explicit quantities, since there is no information
+        // in the WellState to do a decent initialization
+
+        // TODO: Let us assume that the simulator is updated
+
+        // Let us try to do a normal simualtion running, to keep checking the operability status
+        // If the well is not operable during any of the time. It means it does not pass the physical
+        // limit test.
+
+        // create a copy of the well_state to use. If the operability checking is sucessful, we use this one
+        // to replace the original one
+        WellState well_state_copy = well_state;
+
+        // TODO: well state for this well is kind of all zero status
+        // we should be able to provide a better initialization
+        calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger);
+
+        updateWellOperability(ebos_simulator, well_state_copy, deferred_logger);
+
+        if ( !this->isOperable() ) {
+            const std::string msg = " well " + name() + " is not operable during well testing for physical reason";
+            deferred_logger.debug(msg);
+            return;
+        }
+
+        updateWellStateWithTarget(ebos_simulator, well_state_copy, deferred_logger);
+
+        calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger);
+
+        const double dt = ebos_simulator.timeStepSize();
+        const bool converged = this->iterateWellEquations(ebos_simulator, B_avg, dt, well_state_copy, deferred_logger);
+
+        if (!converged) {
+            const std::string msg = " well " + name() + " did not get converged during well testing for physical reason";
+            deferred_logger.debug(msg);
+            return;
+        }
+
+        if (this->isOperable() ) {
+            welltest_state.openWell(name(), WellTestConfig::PHYSICAL );
+            const std::string msg = " well " + name() + " is re-opened through well testing for physical reason";
+            deferred_logger.info(msg);
+            well_state = well_state_copy;
+        } else {
+            const std::string msg = " well " + name() + " is not operable during well testing for physical reason";
+            deferred_logger.debug(msg);
+        }
+    }
+
+
+
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    checkWellOperability(const Simulator& ebos_simulator,
+                         const WellState& well_state,
+                         Opm::DeferredLogger& deferred_logger)
+    {
+
+        const bool checkOperability = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheck);
+        if (!checkOperability) {
+            return;
+        }
+
+        // focusing on PRODUCER for now
+        if (this->isInjector()) {
+            return;
+        }
+
+        if (!this->underPredictionMode() ) {
+            return;
+        }
+
+        if (this->wellIsStopped() && !changed_to_stopped_this_step_) {
+            return;
+        }
+
+        const bool old_well_operable = this->operability_status_.isOperable();
+
+        updateWellOperability(ebos_simulator, well_state, deferred_logger);
+
+        const bool well_operable = this->operability_status_.isOperable();
+
+        if (!well_operable && old_well_operable) {
+            if (well_ecl_.getAutomaticShutIn()) {
+                deferred_logger.info(" well " + name() + " gets SHUT during iteration ");
+            } else {
+                if (!this->wellIsStopped()) {
+                    deferred_logger.info(" well " + name() + " gets STOPPED during iteration ");
+                    this->stopWell();
+                    changed_to_stopped_this_step_ = true;
+                }
+            }
+        } else if (well_operable && !old_well_operable) {
+            deferred_logger.info(" well " + name() + " gets REVIVED during iteration ");
+            this->openWell();
+            changed_to_stopped_this_step_ = false;
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    updateWellOperability(const Simulator& ebos_simulator,
+                          const WellState& well_state,
+                          Opm::DeferredLogger& deferred_logger)
+    {
+        this->operability_status_.reset();
+
+        updateIPR(ebos_simulator, deferred_logger);
+
+        // checking the BHP limit related
+        checkOperabilityUnderBHPLimitProducer(well_state, ebos_simulator, deferred_logger);
+
+        const auto& summaryState = ebos_simulator.vanguard().summaryState();
+
+        // checking whether the well can operate under the THP constraints.
+        if (this->wellHasTHPConstraints(summaryState)) {
+            checkOperabilityUnderTHPLimitProducer(ebos_simulator, well_state, deferred_logger);
+        }
+    }