Support NETBALAN (partially).

NETBALAN can be used to choose balancing frequency (limited: either start
of timestep, or first NUPCOL iterations), balancing tolerance, and maximum
balancing iterations. These correspond to items 1, 2 and 3 of NETBALAN.

opm/simulators/wells/BlackoilWellModel.hpp

@@ -282,7 +282,8 @@ namespace Opm {
             // at the beginning of each time step (Not report step)
             void prepareTimeStep(DeferredLogger& deferred_logger);
             void initPrimaryVariablesEvaluation() const;
-            bool updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls);
+            bool shouldBalanceNetwork(const int reportStepIndex, const int iterationIdx) const;
+            std::pair<bool, double> updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls);
 
             void updateAndCommunicate(const int reportStepIdx,
                                       const int iterationIdx,
@@ -371,7 +372,7 @@ namespace Opm {
                           const double dt);
             void assembleImpl(const int iterationIdx,
                               const double dt,
-                              const int recursion_level,
+                              const std::size_t recursion_level,
                               DeferredLogger& local_deferredLogger);
 
             // called at the end of a time step

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -2092,14 +2092,14 @@ inferLocalShutWells()
     }
 }
 
-bool
+std::pair<bool, double>
 BlackoilWellModelGeneric::
 updateNetworkPressures(const int reportStepIdx)
 {
     // Get the network and return if inactive.
     const auto& network = schedule()[reportStepIdx].network();
     if (!network.active()) {
-        return false;
+        return { false, 0.0 };
     }
     node_pressures_ = WellGroupHelpers::computeNetworkPressures(network,
                                                                 this->wellState(),
@@ -2110,6 +2110,7 @@ updateNetworkPressures(const int reportStepIdx)
 
     // Set the thp limits of wells
     bool active_limit_change = false;
+    double network_imbalance = 0.0;
     for (auto& well : well_container_generic_) {
         // Producers only, since we so far only support the
         // "extended" network model (properties defined by
@@ -2126,11 +2127,12 @@ updateNetworkPressures(const int reportStepIdx)
                 const bool will_switch_to_thp = ws.thp < new_limit;
                 if (thp_is_limit || will_switch_to_thp) {
                     active_limit_change = true;
+                    network_imbalance = std::max(network_imbalance, std::fabs(new_limit - ws.thp));
                 }
             }
         }
     }
-    return active_limit_change;
+    return { active_limit_change, network_imbalance };
 }
 
 void

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -250,7 +250,7 @@ protected:
 
     bool wasDynamicallyShutThisTimeStep(const int well_index) const;
 
-    bool updateNetworkPressures(const int reportStepIdx);
+    std::pair<bool, double> updateNetworkPressures(const int reportStepIdx);
 
     void updateWsolvent(const Group& group,
                         const int reportStepIdx,

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -901,11 +901,11 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     assembleImpl(const int iterationIdx,
                  const double dt,
-                 const int recursion_level,
+                 const std::size_t recursion_level,
                  DeferredLogger& local_deferredLogger)
     {
 
-        const bool network_changed = updateWellControls(local_deferredLogger, /* check group controls */ true); 
+        const auto [network_changed, network_imbalance] = updateWellControls(local_deferredLogger, /* check group controls */ true);
 
         bool alq_updated = false;
         OPM_BEGIN_PARALLEL_TRY_CATCH();
@@ -933,12 +933,10 @@ namespace Opm {
 
         // Maybe do a recursive call to iterate network and well controls.
         if (network_changed) {
-            // Only re-solve network for the first nupcol newton iterations.
-            const int nupcol = schedule()[reportStepIdx].nupcol();
-            if (iterationIdx <= nupcol) {
-                // Limit the number of iterations in the network re-solve.
-                const int max_local_network_iterations = 3;
-                if (recursion_level < max_local_network_iterations) {
+            if (shouldBalanceNetwork(reportStepIdx, iterationIdx)) {
+                const auto& balance = schedule()[reportStepIdx].network_balance();
+                // Iterate if not converged, and number of iterations is not yet max (NETBALAN item 3).
+                if (recursion_level < balance.pressure_max_iter() && network_imbalance > balance.pressure_tolerance()) {
                     assembleImpl(iterationIdx, dt, recursion_level + 1, local_deferredLogger);
                 }
             }
@@ -1469,20 +1467,47 @@ namespace Opm {
     template<typename TypeTag>
     bool
     BlackoilWellModel<TypeTag>::
+    shouldBalanceNetwork(const int reportStepIdx, const int iterationIdx) const
+    {
+        const auto& balance = schedule()[reportStepIdx].network_balance();
+        if (balance.mode() == Network::Balance::CalcMode::TimeStepStart) {
+            return iterationIdx == 0;
+        } else if (balance.mode() == Network::Balance::CalcMode::NUPCOL) {
+            const int nupcol = schedule()[reportStepIdx].nupcol();
+            return iterationIdx < nupcol;
+        } else {
+            // We do not support any other rebalancing modes,
+            // i.e. TimeInterval based rebalancing is not available.
+            // This should be warned about elsewhere, so we choose to
+            // avoid spamming with a warning here.
+            return false;
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    std::pair<bool, double>
+    BlackoilWellModel<TypeTag>::
     updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls)
     {
         // Even if there are no wells active locally, we cannot
         // return as the DeferredLogger uses global communication.
         // For no well active globally we simply return.
-        if( !wellsActive() ) return false;
+        if( !wellsActive() ) return { false, 0.0 };
 
         const int episodeIdx = ebosSimulator_.episodeIndex();
         const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
         const auto& comm = ebosSimulator_.vanguard().grid().comm();
         updateAndCommunicateGroupData(episodeIdx, iterationIdx);
 
-        const bool local_network_changed = updateNetworkPressures(episodeIdx);
+        const auto [local_network_changed, local_network_imbalance]
+            = shouldBalanceNetwork(episodeIdx, iterationIdx) ?
+            updateNetworkPressures(episodeIdx) : std::make_pair(false, 0.0);
         const bool network_changed = comm.sum(local_network_changed);
+        const double network_imbalance = comm.max(local_network_imbalance);
 
         std::set<std::string> switched_wells;
 
@@ -1536,7 +1561,7 @@ namespace Opm {
         const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
         updateWsolvent(fieldGroup, episodeIdx,  this->nupcolWellState());
 
-        return network_changed;
+        return { network_changed, network_imbalance };
     }