adressing several reviewers comments

opm/simulators/wells/BlackoilWellModel.hpp

@@ -413,8 +413,7 @@ template<class Scalar> class WellContributions;
             Scalar gravity_{};
             std::vector<Scalar> depth_{};
             bool alternative_well_rate_init_{};
-            std::optional<double> well_group_thp_calc_;
+            std::map<std::string, Scalar> well_group_thp_calc_;
-
             std::unique_ptr<RateConverterType> rateConverter_{};
             std::map<std::string, std::unique_ptr<AverageRegionalPressureType>> regionalAveragePressureCalculator_{};
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1282,40 +1282,49 @@ namespace Opm {
                     gratTargetFromSales = group_state.grat_sales_target(group.name());
 
                 WGHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff,
-                                                         gratTargetFromSales, nodeName, group_state,
+                                                  gratTargetFromSales, nodeName, group_state,
-                                                         group.has_gpmaint_control(cmode));
+                                                  group.has_gpmaint_control(cmode));
                 const Scalar orig_target = tcalc.groupTarget(ctrl, local_deferredLogger);
 
-                auto mismatch = [&] (auto well_group_thp) {
+                auto mismatch = [&] (auto group_thp) {
-                    Scalar well_group_rate(0.0);
+                    Scalar group_rate(0.0);
                     Scalar rate(0.0);
                     for (auto& well : this->well_container_) {
                         std::string well_name = well->name();
                         auto& ws = well_state.well(well_name);
                         if (group.hasWell(well_name)) {
-                            well->setDynamicThpLimit(well_group_thp);
+                            well->setDynamicThpLimit(group_thp);
                             const Well& well_ecl = this->wells_ecl_[well->indexOfWell()];
                             const auto inj_controls = Well::InjectionControls(0);
                             const auto prod_controls = well_ecl.productionControls(summary_state);
                             well->iterateWellEqWithSwitching(this->simulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger,  false, false);
                             rate = -tcalc.calcModeRateFromRates(ws.surface_rates);
-                            well_group_rate += rate;
+                            group_rate += rate;
                         }
                     }
-                    return (well_group_rate - orig_target)/orig_target;
+                    return (group_rate - orig_target)/orig_target;
                 };
 
-                Scalar min_thp, max_thp;
+                const auto upbranch = network.uptree_branch(nodeName);
-                std::array<Scalar, 2> range_initial;
+                const auto it = this->node_pressures_.find((*upbranch).uptree_node());
+                const Scalar nodal_pressure = it->second;
+                Scalar well_group_thp = nodal_pressure;
+
+                std::optional<Scalar> autochoke_thp;
+                if (auto iter = this->well_group_thp_calc_.find(nodeName); iter != this->well_group_thp_calc_.end()) {
+                    autochoke_thp = this->well_group_thp_calc_.at(nodeName);
+                }
+
                 //Find an initial bracket
-                if (!this->well_group_thp_calc_.has_value()){
+                std::array<Scalar, 2> range_initial;
+                if (!autochoke_thp.has_value()){
+                    Scalar min_thp, max_thp;
                     // Retrieve the terminal pressure of the associated root of the manifold group
                     std::string node_name =  nodeName;
                     while (!network.node(node_name).terminal_pressure().has_value()) {
                         auto branch = network.uptree_branch(node_name).value();
                         node_name = branch.uptree_node();
                     }
-
                     min_thp = network.node(node_name).terminal_pressure().value();
                     WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp);
                     // Narrow down the bracket
@@ -1328,56 +1337,48 @@ namespace Opm {
                     range_initial = {min_thp, max_thp};
                 }
 
-                const auto upbranch = network.uptree_branch(nodeName);
+                if (!autochoke_thp.has_value() || autochoke_thp.value() > nodal_pressure) {
-                const auto it = this->node_pressures_.find((*upbranch).uptree_node());
+                    // The bracket is based on the initial bracket or on a range based on a previous calculated group thp
-                const Scalar nodal_pressure = it->second;
+                    std::array<Scalar, 2> range = autochoke_thp.has_value() ?
-                Scalar well_group_thp = nodal_pressure;
+                        std::array<Scalar, 2>{0.9 * autochoke_thp.value(), 1.1 * autochoke_thp.value()} : range_initial;
-
-                if (!this->well_group_thp_calc_.has_value() || this->well_group_thp_calc_ > nodal_pressure) {
-                    // The bracket is based on the initial bracket or on a range based on a previous calculated common group thp
-                    std::array<Scalar, 2> range;
-                    this->well_group_thp_calc_.has_value() ? 
-                        range =  {0.9 * this->well_group_thp_calc_.value(), 1.1 * this->well_group_thp_calc_.value()} : 
-                        range = range_initial;
-
                     Scalar low, high;
                     std::optional<Scalar> approximate_solution;
                     const Scalar tolerance1 = thp_tolerance;
-                    local_deferredLogger.debug("Using brute force search to bracket the common THP");
+                    local_deferredLogger.debug("Using brute force search to bracket the group THP");
                     const bool finding_bracket = WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger);
 
                     if (approximate_solution.has_value()) {
-                        this->well_group_thp_calc_ = *approximate_solution;
+                        autochoke_thp = *approximate_solution;
-                        local_deferredLogger.debug("Approximate common THP value found: "  + std::to_string(this->well_group_thp_calc_.value()));
+                        local_deferredLogger.debug("Approximate group THP value found: "  + std::to_string(autochoke_thp.value()));
                     } else if (finding_bracket) {
                         const Scalar tolerance2 = thp_tolerance;
                         const int max_iteration_solve = 100;
                         int iteration = 0;
-                        this->well_group_thp_calc_= RegulaFalsiBisection<ThrowOnError>::
+                        autochoke_thp = RegulaFalsiBisection<ThrowOnError>::
                                          solve(mismatch, low, high, max_iteration_solve, tolerance2, iteration);
                         local_deferredLogger.debug(" bracket = [" + std::to_string(low) + ", " + std::to_string(high) + "], " +
                                                    "iteration = " + std::to_string(iteration));
-                        local_deferredLogger.debug("Common THP value = " + std::to_string(this->well_group_thp_calc_.value()));
+                        local_deferredLogger.debug("Group THP value = " + std::to_string(autochoke_thp.value()));
                     } else {
-                        this->well_group_thp_calc_ = {};
+                        autochoke_thp.reset();
-                        local_deferredLogger.debug("Common THP solve failed due to bracketing failure");
+                        local_deferredLogger.debug("Group THP solve failed due to bracketing failure");
                     }
                 }
-                this->well_group_thp_calc_.has_value() ?
+                 if (autochoke_thp.has_value()) {
-                well_group_thp = std::max(this->well_group_thp_calc_.value(), nodal_pressure) : well_group_thp = nodal_pressure;
+                    well_group_thp_calc_[nodeName] = autochoke_thp.value();
+                    // Note: The node pressure of the auto-choke node is set to well_group_thp in computeNetworkPressures()
+                    // and must be larger or equal to the pressure of the uptree node of its branch.
+                    well_group_thp = std::max(autochoke_thp.value(), nodal_pressure);
+                }
 
                 for (auto& well : this->well_container_) {
                     std::string well_name = well->name();
                     if (group.hasWell(well_name)) {
-                        auto& ws = well_state.well(well_name);
+                        well->setDynamicThpLimit(well_group_thp);
-                        if (ws.production_cmode == Opm::WellProducerCMode::THP) {
-                            well->setDynamicThpLimit(well_group_thp);
-                            ws.thp = well_group_thp;
-                        }
                     }
                 }
 
-                // Use the common group THP in computeNetworkPressures
+                // Use the group THP in computeNetworkPressures().
                 group_state.update_well_group_thp(nodeName, well_group_thp);
             }
         }

opm/simulators/wells/GroupState.cpp

@@ -110,7 +110,7 @@ GroupState::update_well_group_thp(const std::string& gname, const double& thp)
 }
 
 template<class Scalar>
-double GroupState<Scalar>::
+Scalar GroupState<Scalar>::
 GroupState::well_group_thp(const std::string& gname) const 
 {
     auto group_iter = this->group_thp.find(gname);

opm/simulators/wells/GroupState.hpp

@@ -50,7 +50,7 @@ public:
     const std::vector<Scalar>& production_rates(const std::string& gname) const;
 
     void update_well_group_thp(const std::string& gname, const double& thp);
-    double well_group_thp(const std::string& gname) const;
+    Scalar well_group_thp(const std::string& gname) const;
 
     bool has_production_reduction_rates(const std::string& gname) const;
     void update_production_reduction_rates(const std::string& gname,
@@ -203,7 +203,7 @@ private:
     std::map<std::string, Scalar> inj_vrep_rate;
     std::map<std::string, Scalar> m_grat_sales_target;
     std::map<std::string, Scalar> m_gpmaint_target;
-    std::map<std::string, double> group_thp;
+    std::map<std::string, Scalar> group_thp;
 
     std::map<std::pair<Phase, std::string>, Group::InjectionCMode> injection_controls;
     WellContainer<GPMaint::State> gpmaint_state;

opm/simulators/wells/WellBhpThpCalculator.cpp

@@ -1052,8 +1052,8 @@ bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
                   Scalar& min_thp, Scalar& max_thp)
 {
     bool bracket_found = false;
-    const int sample_number = 1000;
+    constexpr int sample_number = 1000; 
-    const Scalar interval = 1E5;
+    constexpr Scalar interval = 1E5; 
     Scalar eq_low = eq(min_thp);
     Scalar eq_high = 0.0;
     for (int i = 0; i < sample_number + 1; ++i) {

opm/simulators/wells/WellBhpThpCalculator.hpp

@@ -119,12 +119,6 @@ public:
             const Well& well, 
             const SummaryState& summary_state) const;
 
-  //! \brief Find limits using brute-force solver.
-  static bool bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
-                                const std::array<Scalar, 2>& range,
-                                Scalar& low, Scalar& high,
-                                DeferredLogger& deferred_logger);
-
   //! \brief Find limits using brute-force solver.
   static bool bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
                                 const std::array<Scalar, 2>& range,
@@ -173,6 +167,13 @@ private:
                        std::optional<Scalar>& approximate_solution,
                        DeferredLogger& deferred_logger) const;
 
+ //! \brief Find limits using brute-force solver.
+  static bool bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
+                                const std::array<Scalar, 2>& range,
+                                Scalar& low, Scalar& high,
+                                DeferredLogger& deferred_logger);
+
+
     Scalar findThpFromBhpIteratively(const std::function<Scalar(const Scalar, const Scalar)>& thp_func,
                                      const Scalar bhp,
                                      const Scalar thp_limit,

opm/simulators/wells/WellGroupHelpers.cpp

@@ -937,11 +937,8 @@ computeNetworkPressures(const Network::ExtNetwork& network,
                 } else {
                     // Table number specified as 9999 in the deck, no pressure loss.
                     if (network.node(node).as_choke()){
-                        // Node pressure is set to the common THP of the wells.
+                        // Node pressure is set to the group THP.
-                        // The choke pressure must be non-negative therefore the node pressure of 
+                        node_pressures[node] = group_state.well_group_thp(node);
-                        // the auto-choke node must be larger or equal to the pressure of the uptree node of its branch 
-                        const auto group_thp = group_state.well_group_thp(node);
-                        node_pressures[node] = group_thp >= up_press ? group_thp : up_press;
                     } else {
                             node_pressures[node] = up_press;
                     }

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -487,13 +487,6 @@ WellInterfaceGeneric<Scalar>::getDynamicThpLimit() const
     return dynamic_thp_limit_;
 }
 
-template<class Scalar>
-void WellInterfaceGeneric<Scalar>::
-setDynamicThpLimit(const std::optional<Scalar> thp_limit)
-{
-    dynamic_thp_limit_ = thp_limit;
-}
-
 template<class Scalar>
 void WellInterfaceGeneric<Scalar>::
 updatePerforatedCell(std::vector<bool>& is_cell_perforated)