we dampen the nodal pressure update for each iteration

to improve the network update convergence. The current network update is
somehow explicit.

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -42,6 +42,7 @@
 #include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellTestConfig.hpp>
 #include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
+#include <opm/input/eclipse/Units/Units.hpp>
 
 #include <opm/simulators/utils/DeferredLogger.hpp>
 #include <opm/simulators/wells/BlackoilWellModelConstraints.hpp>
@@ -958,15 +959,18 @@ inferLocalShutWells()
     }
 }
 
-std::pair<bool, double>
+double
 BlackoilWellModelGeneric::
 updateNetworkPressures(const int reportStepIdx)
 {
     // Get the network and return if inactive.
     const auto& network = schedule()[reportStepIdx].network();
     if (!network.active()) {
-        return { false, 0.0 };
+        return 0.0;
     }
+
+    const auto previous_node_pressures = node_pressures_;
+
     node_pressures_ = WellGroupHelpers::computeNetworkPressures(network,
                                                                 this->wellState(),
                                                                 this->groupState(),
@@ -974,10 +978,40 @@ updateNetworkPressures(const int reportStepIdx)
                                                                 schedule(),
                                                                 reportStepIdx);
 
-    // Set the thp limits of wells
+    // here, the network imbalance is the difference between the previous nodal pressure and the new nodal pressure
-    bool active_limit_change = false;
+    double network_imbalance = 0.;
-    double network_imbalance = 0.0;
+
+    if (!previous_node_pressures.empty()) {
+        for (const auto& [name, pressure]: previous_node_pressures) {
+            const auto new_pressure = node_pressures_.at(name);
+            const double change = (new_pressure - pressure);
+            if (std::abs(change) > network_imbalance) {
+                network_imbalance = std::abs(change);
+            }
+            // we dampen the amount of the nodal pressure can change during one iteration
+            // due to the fact our nodal pressure calculation is somewhat explicit
+            // TODO: the following parameters are subject to adjustment for optimization purpose
+            constexpr double upper_update_bound = 5.0 * unit::barsa;
+            constexpr double lower_update_bound = 0.05 * unit::barsa;
+            // relative dampening factor based on update value
+            constexpr double damping_factor = 0.1;
+            const double allowed_change = std::max(std::min(damping_factor * std::abs(change), upper_update_bound),
+                                                   lower_update_bound);
+            if (std::abs(change) > allowed_change) {
+                const double sign = change > 0 ? 1. : -1.;
+                node_pressures_[name] = pressure + sign * allowed_change;
+            }
+        }
+    } else {
+        for (const auto& [name, pressure]: node_pressures_) {
+            if (std::abs(pressure) > network_imbalance) {
+                network_imbalance = std::abs(pressure);
+            }
+        }
+    }
+
     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.
@@ -990,19 +1024,14 @@ updateNetworkPressures(const int reportStepIdx)
                 well->setDynamicThpLimit(new_limit);
                 SingleWellState& ws = this->wellState()[well->indexOfWell()];
                 const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP;
-                const bool will_switch_to_thp = ws.thp < new_limit;
+                // TODO: not sure why the thp is NOT updated properly elsewhere
-                if (thp_is_limit || will_switch_to_thp) {
-                    active_limit_change = true;
-                    network_imbalance = std::max(network_imbalance, std::fabs(new_limit - ws.thp));
-                }
-
                 if (thp_is_limit) {
                     ws.thp = well->getTHPConstraint(summaryState_);
                 }
             }
         }
     }
-    return { active_limit_change, network_imbalance };
+    return network_imbalance;
 }
 
 void

opm/simulators/wells/BlackoilWellModelGeneric.hpp

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

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1573,12 +1573,9 @@ namespace Opm {
         // network related
         bool more_network_update = false;
         if (shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) {
-            const auto [local_network_changed, local_network_imbalance] = updateNetworkPressures(episodeIdx);
+            const auto local_network_imbalance = updateNetworkPressures(episodeIdx);
-            const bool network_changed = comm.sum(local_network_changed);
             const double network_imbalance = comm.max(local_network_imbalance);
-            if (network_changed) {
+            more_network_update = moreNetworkIteration(episodeIdx, network_update_it, network_imbalance);
-                more_network_update = moreNetworkIteration(episodeIdx, network_update_it, network_imbalance);
-            }
         }
 
         bool changed_well_group = false;