Merge pull request #1636 from GitPaean/update_IPR

update well state with THP target based on inflow performance relationship

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -847,7 +847,7 @@ namespace Opm {
         wellhelpers::WellSwitchingLogger logger;
 
         for (const auto& well : well_container_) {
-            well->updateWellControl(well_state_, logger);
+            well->updateWellControl(ebosSimulator_, well_state_, logger);
         }
 
         updateGroupControls();
@@ -938,7 +938,7 @@ namespace Opm {
             well_state_.currentControls()[w] = control;
 
             if (well_state_.effectiveEventsOccurred(w) ) {
-                well->updateWellStateWithTarget(well_state_);
+                well->updateWellStateWithTarget(ebosSimulator_, well_state_);
             }
 
             // there is no new well control change input within a report step,
@@ -1186,7 +1186,7 @@ namespace Opm {
 
             // TODO: we should only do the well is involved in the update group targets
             for (auto& well : well_container_) {
-                well->updateWellStateWithTarget(well_state_);
+                well->updateWellStateWithTarget(ebosSimulator_, well_state_);
                 well->updatePrimaryVariables(well_state_);
             }
         }

opm/autodiff/MultisegmentWell.hpp

@@ -117,9 +117,9 @@ namespace Opm
                                     const double dt,
                                     WellState& well_state) override;
 
-        /// updating the well state based the control mode specified with current
-        // TODO: later will check wheter we need current
-        virtual void updateWellStateWithTarget(WellState& well_state) const override;
+        /// updating the well state based the current control mode
+        virtual void updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
+                                               WellState& well_state) /* const */ override;
 
         /// check whether the well equations get converged for this well
         virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const override;

opm/autodiff/MultisegmentWell_impl.hpp

@@ -254,7 +254,8 @@ namespace Opm
     template <typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    updateWellStateWithTarget(WellState& well_state) const
+    updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
+                              WellState& well_state) /* const */
     {
         // Updating well state bas on well control
         // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE

opm/autodiff/StandardWell.hpp

@@ -140,9 +140,8 @@ namespace Opm
                                     const double dt,
                                     WellState& well_state) override;
 
-        /// updating the well state based the control mode specified with current
-        // TODO: later will check wheter we need current
-        virtual void updateWellStateWithTarget(WellState& well_state) const override;
+        virtual void updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
+                                               WellState& well_state) /* const */ override;
 
         /// check whether the well equations get converged for this well
         virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const override;
@@ -257,6 +256,13 @@ 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_;
+
         const EvalWell& getBhp() const;
 
         EvalWell getQs(const int comp_idx) const;
@@ -354,6 +360,21 @@ 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) const;
+
+        // update WellState based on IPR and associated VFP table
+        void updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator,
+                                             WellState& well_state) const;
+
+        void updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator,
+                                                     WellState& well_state) 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
+        double calculateBHPWithTHPTargetIPR() const;
+
         // relaxation factor considering only one fraction value
         static double relaxationFactorFraction(const double old_value,
                                                const double dx);

opm/autodiff/StandardWell_impl.hpp

@@ -36,6 +36,8 @@ namespace Opm
     , primary_variables_(numWellEq, 0.0)
     , primary_variables_evaluation_(numWellEq) // the number of the primary variables
     , F0_(numWellConservationEq)
+    , ipr_a_(number_of_phases_)
+    , ipr_b_(number_of_phases_)
     {
         assert(num_components_ == numWellConservationEq);
 
@@ -1123,7 +1125,8 @@ namespace Opm
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
-    updateWellStateWithTarget(WellState& well_state) const
+    updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
+                              WellState& well_state) /* const */
     {
         // number of phases
         const int np = number_of_phases_;
@@ -1140,32 +1143,23 @@ namespace Opm
             // TODO: similar to the way below to handle THP
             // we should not something related to thp here when there is thp constraint
             // or when can calculate the THP (table avaiable or requested for output?)
+            // 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: this will be the big task here.
-            // p_bhp = BHP(THP, rates(p_bhp))
-            // more sophiscated techniques is required to obtain the bhp and rates here
-            well_state.thp()[well_index] = target;
-
-            const Opm::PhaseUsage& pu = phaseUsage();
-            std::vector<double> rates(3, 0.0);
-            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-                rates[ Water ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
-            }
-            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
-                 rates[ Oil ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
-            }
-            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                rates[ Gas ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
-            }
-
-            well_state.bhp()[well_index] = calculateBhpFromThp(rates, current);
+            // 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);
             break;
         }
 
         case RESERVOIR_RATE: // intentional fall-through
         case SURFACE_RATE:
+            // TODO: something needs to be done with BHP and THP here
+            // TODO: they should go to a separate function
             // checking the number of the phases under control
             int numPhasesWithTargetsUnderThisControl = 0;
             for (int phase = 0; phase < np; ++phase) {
@@ -1229,6 +1223,196 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    updateIPR(const Simulator& ebos_simulator) 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 (well_type_ == INJECTOR) {
+            return;
+        }
+
+        // initialize all the values to be zero to begin with
+        std::fill(ipr_a_.begin(), ipr_a_.end(), 0.);
+        std::fill(ipr_b_.begin(), ipr_b_.end(), 0.);
+
+        for (int perf = 0; perf < number_of_perforations_; ++perf) {
+            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
+            const double p_r = fs.pressure(FluidSystem::oilPhaseIdx).value();
+
+            // calculating the b for the connection
+            std::vector<double> b_perf(num_components_);
+            for (int 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 = perf_pressure_diffs_[perf];
+            const double pressure_diff = p_r - 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.) {
+                OpmLog::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];
+
+            // 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
+    StandardWell<TypeTag>::
+    updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator,
+                                    WellState& well_state) const
+    {
+        if (well_type_ == PRODUCER) {
+            updateWellStateWithTHPTargetIPRProducer(ebos_simulator,
+                                                    well_state);
+        }
+
+        if (well_type_ == INJECTOR) {
+            well_state.thp()[index_of_well_] = this->getTHPConstraint();
+            // TODO: more work needs to be done for the injectors here, while injectors
+            // have been okay with the current strategy relying on well control equation directly.
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator,
+                                            WellState& well_state) const
+    {
+
+        well_state.thp()[index_of_well_] = this->getTHPConstraint();
+
+        const double bhp = calculateBHPWithTHPTargetIPR();
+
+        well_state.bhp()[index_of_well_] = bhp;
+
+        // TODO: explicit quantities are always tricky for this type of situation
+        updatePrimaryVariables(well_state);
+        initPrimaryVariablesEvaluation();
+
+        std::vector<double> rates;
+        computeWellRatesWithBhp(ebos_simulator, bhp, rates);
+
+        // TODO: double checke the obtained rates
+        // this is another places we might obtain negative rates
+
+        for (size_t p = 0; p < number_of_phases_; ++p) {
+            well_state.wellRates()[number_of_phases_ * index_of_well_ + p] = rates[p];
+        }
+
+        // TODO: there will be something need to be done for the cases not the defaulted 3 phases,
+        // like 2 phases or solvent, polymer, etc. But we are not addressing them with THP control yet.
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    double
+    StandardWell<TypeTag>::
+    calculateBHPWithTHPTargetIPR() const
+    {
+        const double thp_target = this->getTHPConstraint();
+        const double thp_control_index = this->getTHPControlIndex();
+        const  int thp_table_id = well_controls_iget_vfp(well_controls_, thp_control_index);
+        const double alq = well_controls_iget_alq(well_controls_, thp_control_index);
+
+        // not considering injectors for now
+        const double vfp_ref_depth = vfp_properties_->getProd()->getTable(thp_table_id)->getDatumDepth();
+
+        // the density of the top perforation
+        // TODO: make sure this is properly initialized
+        // TODO: with IPR, it should be possible, even this well is a new well, we do not have
+        // TODO: any information of previous rates. However, we are lacking the pressure though.
+        // TODO: but let us not do more work related now to see what will happen
+        const double rho = perf_densities_[0];
+
+        // TODO: call a/the function for dp
+        const double dp = (vfp_ref_depth - ref_depth_) * rho * gravity_;
+
+        const double  bhp_limit = mostStrictBhpFromBhpLimits();
+
+        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;
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::

opm/autodiff/VFPHelpers.hpp

@@ -764,8 +764,93 @@ inline double findTHP(
 
 
 
+// a data type use to do the intersection calculation to get the intial bhp under THP control
+struct RateBhpPair {
+    double rate;
+    double bhp;
+};
 
 
+// looking for a intersection point a line segment and a line, they are both defined with two points
+// it is copied from #include <opm/polymer/Point2D.hpp>, which should be removed since it is only required by the lagacy polymer
+inline bool findIntersection(const std::array<RateBhpPair, 2>& line_segment, const std::array<RateBhpPair, 2>& line, double& bhp) {
+    const double x1 = line_segment[0].rate;
+    const double y1 = line_segment[0].bhp;
+    const double x2 = line_segment[1].rate;
+    const double y2 = line_segment[1].bhp;
+
+    const double x3 = line[0].rate;
+    const double y3 = line[0].bhp;
+    const double x4 = line[1].rate;
+    const double y4 = line[1].bhp;
+
+    const double d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);
+
+    if (d == 0.) {
+        return false;
+    }
+
+    const double x = ((x3 - x4) * (x1 * y2 - y1 * x2) - (x1 - x2) * (x3 * y4 - y3 * x4)) / d;
+    const double y = ((y3 - y4) * (x1 * y2 - y1 * x2) - (y1 - y2) * (x3 * y4 - y3 * x4)) / d;
+
+    if (x >= std::min(x1,x2) && x <= std::max(x1,x2)) {
+        bhp = y;
+        return true;
+    } else {
+        return false;
+    }
+}
+
+// calculating the BHP from thp through the intersection of VFP curves and inflow performance relationship
+inline bool findIntersectionForBhp(const std::vector<RateBhpPair>& ratebhp_samples,
+                                   const std::array<RateBhpPair, 2>& ratebhp_twopoints_ipr,
+                                   double& obtained_bhp)
+{
+    // there possibly two intersection point, then we choose the one corresponding with the bigger rate
+
+    const double bhp1 = ratebhp_twopoints_ipr[0].bhp;
+    const double rate1 = ratebhp_twopoints_ipr[0].rate;
+
+    const double bhp2 = ratebhp_twopoints_ipr[1].bhp;
+    const double rate2 = ratebhp_twopoints_ipr[1].rate;
+
+    assert(rate1 != rate2);
+
+    const double line_slope = (bhp2 - bhp1) / (rate2 - rate1);
+
+    // line equation will be
+    // bhp - bhp1 - line_slope * (flo_rate - flo_rate1) = 0
+    auto flambda = [&](const double flo_rate, const double bhp) {
+        return bhp - bhp1 - line_slope * (flo_rate - rate1);
+    };
+
+    int number_intersection_found = 0;
+    int index_segment = 0; // the intersection segment that intersection happens
+    const size_t num_samples = ratebhp_samples.size();
+    for (size_t i = 0; i < num_samples - 1; ++i) {
+        const double temp1 = flambda(ratebhp_samples[i].rate, ratebhp_samples[i].bhp);
+        const double temp2 = flambda(ratebhp_samples[i+1].rate, ratebhp_samples[i+1].bhp);
+        if (temp1 * temp2 <= 0.) { // intersection happens
+            // in theory there should be maximum two intersection points
+            // while considering the situation == 0. here, we might find more
+            // we always use the last one, which is the one corresponds to the biggest rate,
+            // which we assume is the more stable one
+            ++number_intersection_found;
+            index_segment = i;
+        }
+    }
+
+    if (number_intersection_found == 0) { // there is not intersection point
+        return false;
+    }
+
+    // then we pick the segment from the VFP curve to do the line intersection calculation
+    const std::array<RateBhpPair, 2> line_segment{ ratebhp_samples[index_segment], ratebhp_samples[index_segment + 1] };
+
+    const bool intersection_found = findIntersection(line_segment, ratebhp_twopoints_ipr, obtained_bhp);
+
+    return intersection_found;
+}
 
 
 } // namespace detail

opm/autodiff/VFPProdProperties.cpp

@@ -110,4 +110,147 @@ bool VFPProdProperties::hasTable(const int table_id) const {
 }
 
 
+std::vector<double>
+VFPProdProperties::
+bhpwithflo(const std::vector<double>& flos,
+           const int table_id,
+           const double wfr,
+           const double gfr,
+           const double thp,
+           const double alq,
+           const double dp) const
+{
+    // Get the table
+    const VFPProdTable* table = detail::getTable(m_tables, table_id);
+    const auto thp_i = detail::findInterpData( thp, table->getTHPAxis()); // assume constant
+    const auto wfr_i = detail::findInterpData( wfr, table->getWFRAxis());
+    const auto gfr_i = detail::findInterpData( gfr, table->getGFRAxis());
+    const auto alq_i = detail::findInterpData( alq, table->getALQAxis()); //assume constant
+
+    std::vector<double> bhps(flos.size(), 0.);
+    for (size_t i = 0; i < flos.size(); ++i) {
+        // Value of FLO is negative in OPM for producers, but positive in VFP table
+        const auto flo_i = detail::findInterpData(-flos[i], table->getFloAxis());
+        const detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
+
+        // TODO: this kind of breaks the conventions for the functions here by putting dp within the function
+        bhps[i] = bhp_val.value - dp;
+    }
+
+    return bhps;
+}
+
+
+
+
+
+double
+VFPProdProperties::
+calculateBhpWithTHPTarget(const std::vector<double>& ipr_a,
+                          const std::vector<double>& ipr_b,
+                          const double bhp_limit,
+                          const double thp_table_id,
+                          const double thp_limit,
+                          const double alq,
+                          const double dp) const
+{
+    // For producers, bhp_safe_limit is the highest BHP value that can still produce based on IPR
+    double bhp_safe_limit = 1.e100;
+    for (size_t i = 0; i < ipr_a.size(); ++i) {
+        if (ipr_b[i] == 0.) continue;
+
+        const double bhp = ipr_a[i] / ipr_b[i];
+        if (bhp < bhp_safe_limit) {
+            bhp_safe_limit = bhp;
+        }
+    }
+
+    // Here, we use the middle point between the bhp_limit and bhp_safe_limit to calculate the ratio of the flow
+    // and the middle point serves one of the two points to describe inflow performance relationship line
+    const double bhp_middle = (bhp_limit + bhp_safe_limit) / 2.0;
+
+    // FLO is the rate based on the type specified with the VFP table
+    // The two points correspond to the bhp values of bhp_limit, and the middle of bhp_limit and bhp_safe_limit
+    // for producers, the rates are negative
+    std::vector<double> rates_bhp_limit(ipr_a.size());
+    std::vector<double> rates_bhp_middle(ipr_a.size());
+    for (size_t i = 0; i < rates_bhp_limit.size(); ++i) {
+        rates_bhp_limit[i] = bhp_limit * ipr_b[i] - ipr_a[i];
+        rates_bhp_middle[i] = bhp_middle * ipr_b[i] - ipr_a[i];
+    }
+
+    // TODO: we need to be careful that there is nothings wrong related to the indices here
+    const int Water = BlackoilPhases::Aqua;
+    const int Oil = BlackoilPhases::Liquid;
+    const int Gas = BlackoilPhases::Vapour;
+
+    const VFPProdTable* table = detail::getTable(m_tables, thp_table_id);
+    const double aqua_bhp_limit = rates_bhp_limit[Water];
+    const double liquid_bhp_limit = rates_bhp_limit[Oil];
+    const double vapour_bhp_limit = rates_bhp_limit[Gas];
+    const double flo_bhp_limit = detail::getFlo(aqua_bhp_limit, liquid_bhp_limit, vapour_bhp_limit, table->getFloType() );
+
+    const double aqua_bhp_middle = rates_bhp_middle[Water];
+    const double liquid_bhp_middle = rates_bhp_middle[Oil];
+    const double vapour_bhp_middle = rates_bhp_middle[Gas];
+    const double flo_bhp_middle = detail::getFlo(aqua_bhp_middle, liquid_bhp_middle, vapour_bhp_middle, table->getFloType() );
+
+    // we use the ratios based on the middle value of bhp_limit and bhp_safe_limit
+    const double wfr = detail::getWFR(aqua_bhp_middle, liquid_bhp_middle, vapour_bhp_middle, table->getWFRType());
+    const double gfr = detail::getGFR(aqua_bhp_middle, liquid_bhp_middle, vapour_bhp_middle, table->getGFRType());
+
+    // we get the flo sampling points from the table,
+    // then extend it with zero and rate under bhp_limit for extrapolation
+    std::vector<double> flo_samples = table->getFloAxis();
+
+    if (flo_samples[0] > 0.) {
+        flo_samples.insert(flo_samples.begin(), 0.);
+    }
+
+    if (flo_samples.back() < std::abs(flo_bhp_limit)) {
+        flo_samples.push_back(std::abs(flo_bhp_limit));
+    }
+
+    // kind of unncessarily following the tradation that producers should have negative rates
+    // the key is here that it should be consistent with the function bhpwithflo
+    for (double& value : flo_samples) {
+        value = -value;
+    }
+
+    // get the bhp sampling values based on the flo sample values
+    const std::vector<double> bhp_flo_samples = bhpwithflo(flo_samples, thp_table_id, wfr, gfr, thp_limit, alq, dp);
+
+    std::vector<detail::RateBhpPair> ratebhp_samples;
+    for (size_t i = 0; i < flo_samples.size(); ++i) {
+        ratebhp_samples.push_back( detail::RateBhpPair{flo_samples[i], bhp_flo_samples[i]} );
+    }
+
+    const std::array<detail::RateBhpPair, 2> ratebhp_twopoints_ipr {detail::RateBhpPair{flo_bhp_middle, bhp_middle},
+                                                                    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);
+
+    // \Note: assuming not that negative BHP does not make sense
+    if (obtain_solution_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 "
+                                   " cause problems later ");
+            OpmLog::info("TOO_HIGH_BHP_FOUND_THP_TARGET", msg);
+
+            const std::string debug_msg = " obtain_bhp " + std::to_string(obtain_bhp)
+                                        + " bhp_safe_limit " + std::to_string(bhp_safe_limit)
+                                        + " thp limit " + std::to_string(thp_limit);
+            OpmLog::debug(debug_msg);
+        }
+        return obtain_bhp;
+    } else {
+        OpmLog::warning("NO_BHP_FOUND_THP_TARGET", " we could not find a bhp value with thp target.");
+        return -100.;
+    }
+}
+
+
+
 }

opm/autodiff/VFPProdProperties.hpp

@@ -170,7 +170,30 @@ public:
         return m_tables.empty();
     }
 
+
+    /**
+     * Calculate the Bhp value from the THP target/constraint value
+     * based on inflow performance relationship and VFP curves
+     */
+     double
+     calculateBhpWithTHPTarget(const std::vector<double>& ipr_a,
+                               const std::vector<double>& ipr_b,
+                               const double bhp_limit,
+                               const double thp_table_id,
+                               const double thp_limit,
+                               const double alq,
+                               const double dp) const;
+
 protected:
+    // calculate a group bhp values with a group of flo rate values
+    std::vector<double> bhpwithflo(const std::vector<double>& flos,
+                                   const int table_id,
+                                   const double wfr,
+                                   const double gfr,
+                                   const double thp,
+                                   const double alq,
+                                   const double dp) const;
+
     // Map which connects the table number with the table itself
     std::map<int, const VFPProdTable*> m_tables;
 };

opm/autodiff/WellInterface.hpp

@@ -172,10 +172,12 @@ namespace Opm
                                            const WellState& well_state,
                                            std::vector<double>& well_potentials) = 0;
 
-        virtual void updateWellStateWithTarget(WellState& well_state) const = 0;
+        virtual void updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
+                                               WellState& well_state) /* const */ = 0;
 
-        void updateWellControl(WellState& well_state,
-                               wellhelpers::WellSwitchingLogger& logger) const;
+        void updateWellControl(/* const */ Simulator& ebos_simulator,
+                               WellState& well_state,
+                               wellhelpers::WellSwitchingLogger& logger) /* const */;
 
         virtual void updatePrimaryVariables(const WellState& well_state) const = 0;
 
@@ -312,8 +314,14 @@ namespace Opm
         bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
                                  const WellState& well_state) const;
 
+        bool underPredictionMode() const;
+
         bool wellHasTHPConstraints() const;
 
+        double getTHPConstraint() const;
+
+        int getTHPControlIndex() const;
+
         // Component fractions for each phase for the well
         const std::vector<double>& compFrac() const;
 

opm/autodiff/WellInterface_impl.hpp

@@ -364,18 +364,48 @@ namespace Opm
 
 
 
+
     template<typename TypeTag>
     bool
     WellInterface<TypeTag>::
     wellHasTHPConstraints() const
+    {
+        return getTHPControlIndex() >= 0;
+    }
+
+
+
+
+    template<typename TypeTag>
+    double
+    WellInterface<TypeTag>::
+    getTHPConstraint() const
+    {
+        const int thp_control_index = getTHPControlIndex();
+
+        if (thp_control_index < 0) {
+            OPM_THROW(std::runtime_error, " there is no THP constraint/limit for well " << name()
+                                          << ", while we are requesting it ");
+        }
+
+        return well_controls_iget_target(well_controls_, thp_control_index);
+    }
+
+
+
+
+    template<typename TypeTag>
+    int
+    WellInterface<TypeTag>::
+    getTHPControlIndex() const
     {
         const int nwc = well_controls_get_num(well_controls_);
         for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
             if (well_controls_iget_type(well_controls_, ctrl_index) == THP) {
-                return true;
+                return ctrl_index;
             }
         }
-        return false;
+        return -1;
     }
 
 
@@ -385,8 +415,9 @@ namespace Opm
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::
-    updateWellControl(WellState& well_state,
-                      wellhelpers::WellSwitchingLogger& logger) const
+    updateWellControl(/* const */ Simulator& ebos_simulator,
+                      WellState& well_state,
+                      wellhelpers::WellSwitchingLogger& logger) /* const */
     {
         const int np = number_of_phases_;
         const int w = index_of_well_;
@@ -437,7 +468,7 @@ namespace Opm
         }
 
         if (updated_control_index != old_control_index) { //  || well_collection_->groupControlActive()) {
-            updateWellStateWithTarget(well_state);
+            updateWellStateWithTarget(ebos_simulator, well_state);
             updatePrimaryVariables(well_state);
         }
     }
@@ -446,6 +477,31 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    bool
+    WellInterface<TypeTag>::
+    underPredictionMode() const
+    {
+        bool under_prediction_mode = false;
+
+        switch( well_type_ ) {
+        case PRODUCER:
+            under_prediction_mode = well_ecl_->getProductionProperties(current_step_).predictionMode;
+            break;
+        case INJECTOR:
+            under_prediction_mode = well_ecl_->getInjectionProperties(current_step_).predictionMode;
+            break;
+        default:
+            OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name());
+        }
+
+        return under_prediction_mode;
+    }
+
+
+
+
+
     template<typename TypeTag>
     bool
     WellInterface<TypeTag>::
@@ -1078,7 +1134,7 @@ namespace Opm
             solveEqAndUpdateWellState(well_state);
 
             wellhelpers::WellSwitchingLogger logger;
-            updateWellControl(well_state, logger);
+            updateWellControl(ebosSimulator, well_state, logger);
             initPrimaryVariablesEvaluation();
         } while (it < max_iter);