adding function calculateBhpWithTHPTarget to VFPProdProperties

it calculate bhp value based on THP target/limit, VFP curves and
inflow-performance relationship

opm/autodiff/StandardWell_impl.hpp

@@ -1315,8 +1315,9 @@ namespace Opm
             }
 
             for (int p = 0; p < number_of_phases_; ++p) {
-                ipr_a_[p] += ipr_a_perf[p];
+                // TODO: double check the indices here
-                ipr_b_[p] += ipr_b_perf[p];
+                ipr_a_[ebosCompIdxToFlowCompIdx(p)] += ipr_a_perf[p];
+                ipr_b_[ebosCompIdxToFlowCompIdx(p)] += ipr_b_perf[p];
             }
         }
     }

opm/autodiff/VFPHelpers.hpp

@@ -808,7 +808,7 @@ inline bool findIntersectionForBhp(const std::vector<double>&rate_samples,
                                    const double flo_rate2,
                                    const double bhp1,
                                    const double bhp2,
-                                   double& bhp)
+                                   double& obtained_bhp)
 {
     // there possibly two intersection point, then we choose the bigger one
     // we choose the bigger one, then it will be the later one in the rate_samples
@@ -851,10 +851,10 @@ inline bool findIntersectionForBhp(const std::vector<double>&rate_samples,
     const std::array<DataPoint, 2> line { DataPoint{flo_rate1, bhp1},
                                           DataPoint{flo_rate2, bhp2} };
 
-    const bool inter_section_found = findIntersection(line_segment, line, bhp);
+    const bool inter_section_found = findIntersection(line_segment, line, obtained_bhp);
 
     if (inter_section_found) {
-        std::cout << " found a bhp is " << bhp << std::endl;
+        std::cout << " found a bhp is " << obtained_bhp << std::endl;
         return true;
     } else {
         std::cout << " did not find the intersection point " << std::endl;

opm/autodiff/VFPProdProperties.cpp

@@ -110,4 +110,136 @@ 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);
+
+    double obtain_bhp = 0.;
+    const bool obtain_solution_with_thp_limit = detail::findIntersectionForBhp(flo_samples, bhp_flo_samples,
+                                                        flo_bhp_middle, flo_bhp_limit, bhp_middle, bhp_limit, 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) {
+            std::cout << " Look like we are getting a too high BHP value from the THP constraint "
+                      << " which might cause problems later " << std::endl;
+
+            std::cout << " obtain_bhp " << obtain_bhp << " bhp_safe_limit " << bhp_safe_limit << std::endl;
+        }
+        return obtain_bhp;
+    } else {
+        std::cout << " COULD NOT find an Intersection point " << std::endl;
+        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;
 };