allow individual well constraints

before rebasing

moved common thp calculation to updateWellControls

Small fix

clean up and improvements according reviewer comments

clean up and improvements according reviewer comments, part 2

changed assessing safe THP range

rebasing fixes

removed unused argument

rebasing

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1546,6 +1546,7 @@ updateNetworkPressures(const int reportStepIdx, const Scalar damping_factor)
     }
 
     for (auto& well : well_container_generic_) {
+
         // Producers only, since we so far only support the
         // "extended" network model (properties defined by
         // BRANPROP and NODEPROP) which only applies to producers.

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -55,6 +55,7 @@
 #include <opm/simulators/linalg/gpubridge/WellContributions.hpp>
 #endif
 
+
 #if HAVE_MPI
 #include <opm/simulators/utils/MPISerializer.hpp>
 #endif
@@ -1234,9 +1235,6 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     updateWellControlsAndNetwork(const bool mandatory_network_balance, const double dt, DeferredLogger& local_deferredLogger)
     {
-        // PJPE: calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES)
-        computeWellGroupThp(local_deferredLogger);
-
         // not necessarily that we always need to update once of the network solutions
         bool do_network_update = true;
         bool well_group_control_changed = false;
@@ -1331,6 +1329,7 @@ namespace Opm {
         const auto& balance = this->schedule()[reportStepIdx].network_balance();
         const Scalar thp_tolerance = balance.thp_tolerance();
 
+
         if (!network.active()) {
             return;
         }
@@ -1346,6 +1345,7 @@ namespace Opm {
                 const auto ctrl = group.productionControls(summary_state);
                 const auto cmode = ctrl.cmode;
                 const auto pu = this->phase_usage_;
+
                 //TODO: Auto choke combined with RESV control is not supported
                 std::vector<Scalar> resv_coeff(pu.num_phases, 1.0);
                 Scalar gratTargetFromSales = 0.0;
@@ -1376,6 +1376,31 @@ namespace Opm {
                     return (group_rate - orig_target)/orig_target;
                 };
 
+                double min_thp, max_thp;
+                std::array<double, 2> range_initial;
+                //Find an initial bracket
+                if (!this->well_group_thp_calc_.has_value()){
+                    // 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();
+                    std::optional<double> approximate_solution0;
+                    WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp, local_deferredLogger);
+
+                     // Narrow down the bracket
+                    double low1, high1;
+                    std::array<double, 2> range = {0.9*min_thp, 1.1*max_thp};
+                    std::optional<double> appr_sol;
+                    WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low1, high1, appr_sol, 0.0, local_deferredLogger);
+                    min_thp = low1;
+                    max_thp = high1;
+                    range_initial = {min_thp, max_thp};
+                }
+
                 const auto upbranch = network.uptree_branch(nodeName);
                 const auto it = this->node_pressures_.find((*upbranch).uptree_node());
                 const Scalar nodal_pressure = it->second;
@@ -1446,6 +1471,7 @@ namespace Opm {
                 for (auto& well : this->well_container_) {
                     std::string well_name = well->name();
                     if (group.hasWell(well_name)) {
+
                         well->setDynamicThpLimit(well_group_thp);
                     }
                 }
@@ -2686,10 +2712,10 @@ namespace Opm {
                     const std::string msg = "Compute initial well solution for " + well->name() + " initially failed. Continue with the previous rates";
                     deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
                 }
-                // If we're using local well solves that include control switches, they also update
-                // operability, so reset before main iterations begin
-                well->resetWellOperability();
             }
+            // If we're using local well solves that include control switches, they also update
+            // operability, so reset before main iterations begin
+            well->resetWellOperability();
         }
         updatePrimaryVariables(deferred_logger);
 

opm/simulators/wells/GroupState.cpp

@@ -105,8 +105,9 @@ GroupState<Scalar>::production_rates(const std::string& gname) const
 }
 
 //-------------------------------------------------------------------------
+
 template<class Scalar>
-bool GroupState<Scalar>::
+void GroupState<Scalar>::
 GroupState::update_well_group_thp(const std::string& gname, const double& thp) 
 {
     this->group_thp[gname] = thp;
@@ -123,6 +124,8 @@ GroupState::well_group_thp(const std::string& gname) const
     return group_iter->second;
 }
 
+//-------------------------------------------------------------------------
+
 template<class Scalar>
 void GroupState<Scalar>::
 GroupState::update_well_group_thp(const std::string& gname, const double& thp) 

opm/simulators/wells/WellBhpThpCalculator.hpp

@@ -130,6 +130,24 @@ public:
   static bool bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
                                 Scalar& min_thp, Scalar& max_thp);
 
+  //! \brief Find limits using brute-force solver.
+  static bool bruteForceBracket(const std::function<double(const double)>& eq,
+                                const std::array<double, 2>& range,
+                                double& low, double& high,
+                                DeferredLogger& deferred_logger);
+
+  //! \brief Find limits using brute-force solver.
+  static bool bruteForceBracketCommonTHP(const std::function<double(const double)>& eq,
+                                const std::array<double, 2>& range,
+                                double& low, double& high,
+                                std::optional<double>& approximate_solution,
+                                const double& limit,
+                                DeferredLogger& deferred_logger);
+
+  //! \brief Find limits using brute-force solver.
+  static bool bruteForceBracketCommonTHP(const std::function<double(const double)>& eq,
+                                double& min_thp, double& max_thp);
+
 private:
     //! \brief Compute BHP from THP limit for an injector - implementation.
     template<class ErrorPolicy>

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -518,6 +518,13 @@ WellInterfaceGeneric<Scalar>::getDynamicThpLimit() const
     return dynamic_thp_limit_;
 }
 
+template<class Scalar>
+void WellInterfaceGeneric::
+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)

opm/simulators/wells/WellInterface_impl.hpp

@@ -1480,11 +1480,11 @@ namespace Opm
                     }
                     ws.trivial_target = false;
                 } else {
-                        ws.trivial_target = true;
+                    ws.trivial_target = true;
                 }
                     
                 break;
-            }   
+            }
             case Well::ProducerCMode::CMODE_UNDEFINED:
             case Well::ProducerCMode::NONE:
             {