Actually respect NUPCOL (not NUPCOL+1) and allow network to converge also with gas lift optimization.

opm/simulators/wells/BlackoilWellModel.hpp

@@ -74,6 +74,7 @@
 #include <map>
 #include <memory>
 #include <string>
+#include <tuple>
 #include <vector>
 
 namespace Opm {
@@ -267,7 +268,7 @@ template<class Scalar> class WellContributions;
             // at the beginning of each time step (Not report step)
             void prepareTimeStep(DeferredLogger& deferred_logger);
 
-            std::pair<bool, bool>
+            std::tuple<bool, bool, Scalar>
             updateWellControls(const bool mandatory_network_balance, DeferredLogger& deferred_logger, const bool relax_network_tolerance = false);
 
             void updateAndCommunicate(const int reportStepIdx,
@@ -416,8 +417,9 @@ template<class Scalar> class WellContributions;
             // the function handles one iteration of updating well controls and network pressures.
             // it is possible to decouple the update of well controls and network pressures further.
             // the returned two booleans are {continue_due_to_network, well_group_control_changed}, respectively
-            std::pair<bool, bool> updateWellControlsAndNetworkIteration(const bool mandatory_network_balance,
+            std::tuple<bool, bool, Scalar> updateWellControlsAndNetworkIteration(const bool mandatory_network_balance,
                                                                         const bool relax_network_tolerance,
+                                                                        const bool optimize_gas_lift,
                                                                         const double dt,
                                                                         DeferredLogger& local_deferredLogger);
 

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1270,7 +1270,7 @@ updateAndCommunicateGroupData(const int reportStepIdx,
     // before we copy to well_state_nupcol_.
     this->wellState().updateGlobalIsGrup(comm_);
 
-    if (iterationIdx <= nupcol) {
+    if (iterationIdx < nupcol) {
         OPM_TIMEBLOCK(updateNupcol);
         this->updateNupcolWGState();
     } else {

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1192,6 +1192,7 @@ namespace Opm {
         // not necessarily that we always need to update once of the network solutions
         bool do_network_update = true;
         bool well_group_control_changed = false;
+        Scalar network_imbalance = 0.0;
         // after certain number of the iterations, we use relaxed tolerance for the network update
         const std::size_t iteration_to_relax = param_.network_max_strict_outer_iterations_;
         // after certain number of the iterations, we terminate
@@ -1203,12 +1204,16 @@ namespace Opm {
                 const int episodeIdx = simulator_.episodeIndex();
                 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
                 if (this->shouldBalanceNetwork(episodeIdx, iterationIdx + 1)) {
-                    local_deferredLogger.debug("maximum of " + std::to_string(max_iteration) + " network iterations has been used, we stop the update \n "
-                                                " and try again after the next newton iteration.");
+                    const std::string msg = fmt::format("Maximum of {:d} network iterations has been used and we stop the update, \n"
+                        "and try again after the next Newton iteration (imbalance = {:.2e} bar, ctrl_change = {})",
+                        max_iteration, network_imbalance*1.0e-5, well_group_control_changed);
+                    local_deferredLogger.debug(msg);
                 } else {
                     if (this->terminal_output_) {
-                        local_deferredLogger.info("maximum of " + std::to_string(max_iteration) + " network iterations has been used and we stop the update. "
-                                                  "The simulator will continue with unconverged network results.");
+                        const std::string msg = fmt::format("Maximum of {:d} network iterations has been used and we stop the update. \n"
+                            "The simulator will continue with unconverged network results (imbalance = {:.2e} bar, ctrl_change = {})",
+                            max_iteration, network_imbalance*1.0e-5, well_group_control_changed);
+                        local_deferredLogger.info(msg);
                     }
                 }
                 break;
@@ -1217,8 +1222,10 @@ namespace Opm {
                 local_deferredLogger.debug("We begin using relaxed tolerance for network update now after " + std::to_string(iteration_to_relax) + " iterations ");
             }
             const bool relax_network_balance = network_update_iteration >= iteration_to_relax;
-            std::tie(do_network_update, well_group_control_changed) =
-                    updateWellControlsAndNetworkIteration(mandatory_network_balance, relax_network_balance, dt,local_deferredLogger);
+            // Never optimize gas lift in last iteration, to allow network convergence (unless max_iter < 2)
+            const bool optimize_gas_lift = ( (network_update_iteration + 1) < std::max(max_iteration, static_cast<std::size_t>(2)) );
+            std::tie(well_group_control_changed, do_network_update, network_imbalance) =
+                    updateWellControlsAndNetworkIteration(mandatory_network_balance, relax_network_balance, optimize_gas_lift, dt,local_deferredLogger);
             ++network_update_iteration;
         }
         return well_group_control_changed;
@@ -1228,15 +1235,16 @@ namespace Opm {
 
 
     template<typename TypeTag>
-    std::pair<bool, bool>
+    std::tuple<bool, bool, typename BlackoilWellModel<TypeTag>::Scalar>
     BlackoilWellModel<TypeTag>::
     updateWellControlsAndNetworkIteration(const bool mandatory_network_balance,
                                           const bool relax_network_tolerance,
+                                          const bool optimize_gas_lift,
                                           const double dt,
                                           DeferredLogger& local_deferredLogger)
     {
         OPM_TIMEFUNCTION();
-        auto [well_group_control_changed, more_inner_network_update] =
+        auto [well_group_control_changed, more_inner_network_update, network_imbalance] =
                 updateWellControls(mandatory_network_balance,
                                    local_deferredLogger,
                                    relax_network_tolerance);
@@ -1244,7 +1252,7 @@ namespace Opm {
         bool alq_updated = false;
         OPM_BEGIN_PARALLEL_TRY_CATCH();
         {
-            alq_updated = gaslift_.maybeDoGasLiftOptimize(simulator_,
+            if (optimize_gas_lift) alq_updated = gaslift_.maybeDoGasLiftOptimize(simulator_,
                                                           well_container_,
                                                           this->wellState(),
                                                           this->groupState(),
@@ -1283,7 +1291,7 @@ namespace Opm {
         const int iterationIdx = simulator_.model().newtonMethod().numIterations();
         const bool more_network_update = this->shouldBalanceNetwork(reportStepIdx, iterationIdx) &&
                     (more_inner_network_update || well_group_control_changed || alq_updated);
-        return {more_network_update, well_group_control_changed};
+        return {well_group_control_changed, more_network_update, network_imbalance};
     }
 
     // This function is to be used for well groups in an extended network that act as a subsea manifold
@@ -1743,7 +1751,7 @@ namespace Opm {
 
 
     template<typename TypeTag>
-    std::pair<bool, bool>
+    std::tuple<bool, bool, typename BlackoilWellModel<TypeTag>::Scalar>
     BlackoilWellModel<TypeTag>::
     updateWellControls(const bool mandatory_network_balance,
                        DeferredLogger& deferred_logger,
@@ -1753,7 +1761,7 @@ namespace Opm {
         const int episodeIdx = simulator_.episodeIndex();
         const auto& network = this->schedule()[episodeIdx].network();
         if (!this->wellsActive() && !network.active()) {
-            return {false, false};
+            return {false, false, 0.0};
         }
 
         const int iterationIdx = simulator_.model().newtonMethod().numIterations();
@@ -1761,6 +1769,7 @@ namespace Opm {
         this->updateAndCommunicateGroupData(episodeIdx, iterationIdx, param_.nupcol_group_rate_tolerance_, deferred_logger);
 
         // network related
+        Scalar network_imbalance = 0.0;
         bool more_network_update = false;
         if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) {
             OPM_TIMEBLOCK(BalanceNetowork);
@@ -1773,7 +1782,7 @@ namespace Opm {
             bool more_network_sub_update = false;
             for (int i = 0; i < max_number_of_sub_iterations; i++) {
                 const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx, network_pressure_update_damping_factor, network_max_pressure_update);
-                const Scalar network_imbalance = comm.max(local_network_imbalance);
+                network_imbalance = comm.max(local_network_imbalance);
                 const auto& balance = this->schedule()[episodeIdx].network_balance();
                 constexpr Scalar relaxation_factor = 10.0;
                 const Scalar tolerance = relax_network_tolerance ? relaxation_factor * balance.pressure_tolerance() : balance.pressure_tolerance();
@@ -1854,7 +1863,7 @@ namespace Opm {
         // update wsolvent fraction for REIN wells
         this->updateWsolvent(fieldGroup, episodeIdx,  this->nupcolWellState());
 
-        return { changed_well_group, more_network_update };
+        return { changed_well_group, more_network_update, network_imbalance };
     }
 
     template<typename TypeTag>
@@ -1903,7 +1912,7 @@ namespace Opm {
         bool changed = false;
         // restrict the number of group switches but only after nupcol iterations.
         const int nupcol = this->schedule()[reportStepIdx].nupcol();
-        const int max_number_of_group_switches = iterationIdx <= nupcol ? 9999 : param_.max_number_of_group_switches_;
+        const int max_number_of_group_switches = iterationIdx < nupcol ? 9999 : param_.max_number_of_group_switches_;
         bool changed_hc = this->checkGroupHigherConstraints( group, deferred_logger, reportStepIdx, max_number_of_group_switches);
         if (changed_hc) {
             changed = true;

opm/simulators/wells/GasLiftGroupInfo.cpp

@@ -418,10 +418,10 @@ checkNewtonIterationIdxOk_(const std::string& well_name)
                 "LIFTOPT item4 == YES, it = {}, nupcol = {} -->  GLIFT optimize = {}",
                 this->iteration_idx_,
                 nupcol,
-                ((this->iteration_idx_ <= nupcol) ? "TRUE" : "FALSE"));
+                ((this->iteration_idx_ < nupcol) ? "TRUE" : "FALSE"));
             displayDebugMessage_(msg, well_name);
         }
-        return this->iteration_idx_ <= nupcol;
+        return this->iteration_idx_ < nupcol;
     }
     else {
         if (this->debug) {