Merge pull request #4409 from GitPaean/update_well_state_with_thp_target2

update well state when switching to THP control for producers

opm/simulators/wells/MultisegmentWell.hpp

@@ -258,6 +258,10 @@ namespace Opm
                                  const Simulator& ebos_simulator,
                                  DeferredLogger& deferred_logger) const;
 
+        bool updateWellStateWithTHPTargetProd(const Simulator& ebos_simulator,
+                                              WellState& well_state,
+                                              DeferredLogger& deferred_logger) const;
+
         virtual double getRefDensity() const override;
 
         virtual bool iterateWellEqWithControl(const Simulator& ebosSimulator,

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -2020,4 +2020,32 @@ namespace Opm
         connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
     }
 
+
+
+
+    template<typename TypeTag>
+    bool
+    MultisegmentWell<TypeTag>::
+    updateWellStateWithTHPTargetProd(const Simulator& ebos_simulator,
+                                     WellState& well_state,
+                                     DeferredLogger& deferred_logger) const
+    {
+        const auto& summary_state = ebos_simulator.vanguard().summaryState();
+
+        auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
+            ebos_simulator, summary_state, this->getALQ(well_state), deferred_logger);
+        if (bhp_at_thp_limit) {
+            std::vector<double> rates(this->number_of_phases_, 0.0);
+            computeWellRatesWithBhpIterations(ebos_simulator, *bhp_at_thp_limit, rates, deferred_logger);
+            auto& ws = well_state.well(this->name());
+            ws.surface_rates = rates;
+            ws.bhp = *bhp_at_thp_limit;
+            ws.thp = this->getTHPConstraint(summary_state);
+            return true;
+        } else {
+            return false;
+        }
+    }
+
+
 } // namespace Opm

opm/simulators/wells/StandardWell.hpp

@@ -334,6 +334,9 @@ namespace Opm
             DeferredLogger& deferred_logger,
             const WellState &well_state) const;
 
+        bool updateWellStateWithTHPTargetProd(const Simulator& ebos_simulator,
+                                              WellState& well_state,
+                                              DeferredLogger& deferred_logger) const;
 
         virtual double getRefDensity() const override;
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -1783,6 +1783,34 @@ namespace Opm
         return potentials;
     }
 
+
+
+    template<typename TypeTag>
+    bool
+    StandardWell<TypeTag>::
+    updateWellStateWithTHPTargetProd(const Simulator& ebos_simulator,
+                                     WellState& well_state,
+                                     DeferredLogger& deferred_logger) const
+    {
+        const auto& summary_state = ebos_simulator.vanguard().summaryState();
+
+        auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
+            ebos_simulator, summary_state, this->getALQ(well_state), deferred_logger);
+        if (bhp_at_thp_limit) {
+            std::vector<double> rates(this->number_of_phases_, 0.0);
+            computeWellRatesWithBhp(ebos_simulator, *bhp_at_thp_limit, rates, deferred_logger);
+            auto& ws = well_state.well(this->name());
+            ws.surface_rates = rates;
+            ws.bhp = *bhp_at_thp_limit;
+            ws.thp = this->getTHPConstraint(summary_state);
+            return true;
+        } else {
+            return false;
+        }
+    }
+
+
+
     template<typename TypeTag>
     double
     StandardWell<TypeTag>::
@@ -2363,9 +2391,14 @@ namespace Opm
                                                                                alq_value,
                                                                                this->getTHPConstraint(summary_state),
                                                                                deferred_logger);
-        auto v = frates(*bhpAtLimit);
+
-        if (bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }))
+        if (bhpAtLimit) {
-            return bhpAtLimit;
+            auto v = frates(*bhpAtLimit);
+            if (std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }) ) {
+                return bhpAtLimit;
+            }
+        }
+
 
         auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
             // Solver the well iterations to see if we are
@@ -2384,9 +2417,15 @@ namespace Opm
                                                                           alq_value,
                                                                           this->getTHPConstraint(summary_state),
                                                                           deferred_logger);
-        v = frates(*bhpAtLimit);
+
-        if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }))
+
-            return bhpAtLimit;
+        if (bhpAtLimit) {
+            // should we use fratesIter here since fratesIter is used in computeBhpAtThpLimitProd above?
+            auto v = frates(*bhpAtLimit);
+            if (std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }) ) {
+                return bhpAtLimit;
+            }
+        }
 
         // we still don't get a valied solution.
         return std::nullopt;

opm/simulators/wells/WellBhpThpCalculator.cpp

@@ -738,7 +738,7 @@ bruteForceBracket(const std::function<double(const double)>& eq,
     bool finding_bracket = false;
     low = range[0];
     high = range[1];
-    const int sample_number = 100;
+    const int sample_number = 200;
     const double interval = (high - low) / sample_number;
     double eq_low = eq(low);
     double eq_high;

opm/simulators/wells/WellInterface.hpp

@@ -201,6 +201,10 @@ public:
                                            WellState& well_state,
                                            DeferredLogger& deferred_logger) const;
 
+    virtual bool updateWellStateWithTHPTargetProd(const Simulator& ebos_simulator,
+                                                  WellState& well_state,
+                                                  DeferredLogger& deferred_logger) const = 0;
+
     enum class IndividualOrGroup { Individual, Group, Both };
     bool updateWellControl(const Simulator& ebos_simulator,
                            const IndividualOrGroup iog,

opm/simulators/wells/WellInterface_impl.hpp

@@ -965,24 +965,26 @@ namespace Opm
             }
             case Well::ProducerCMode::THP:
             {
-                auto rates = ws.surface_rates;
+                const bool update_success = updateWellStateWithTHPTargetProd(ebos_simulator, well_state, deferred_logger);
-                this->adaptRatesForVFP(rates);
-                double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
-                                                                             rates,
-                                                                             well,
-                                                                             summaryState,
-                                                                             this->getRefDensity(),
-                                                                             deferred_logger);
-                ws.bhp = bhp;
-                ws.thp = this->getTHPConstraint(summaryState);
 
-                // if the total rates are negative or zero
+                if (!update_success) {
-                // we try to provide a better intial well rate
+                    // the following is the original way of initializing well state with THP constraint
-                // using the well potentials
+                    // keeping it for robust reason in case that it fails to get a bhp value with THP constraint
-                double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
+                    // more sophisticated design might be needed in the future
-                if (total_rate <= 0.0){
+                    auto rates = ws.surface_rates;
-                    for (int p = 0; p<np; ++p) {
+                    this->adaptRatesForVFP(rates);
-                        ws.surface_rates[p] = -ws.well_potentials[p];
+                    const double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(
+                        well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
+                    ws.bhp = bhp;
+                    ws.thp = this->getTHPConstraint(summaryState);
+                    // if the total rates are negative or zero
+                    // we try to provide a better initial well rate
+                    // using the well potentials
+                    const double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
+                    if (total_rate <= 0.0) {
+                        for (int p = 0; p < this->number_of_phases_; ++p) {
+                            ws.surface_rates[p] = -ws.well_potentials[p];
+                        }
                     }
                 }
                 break;