Merge pull request #3954 from atgeirr/iterate-network

Add iterations to the network model

opm/simulators/utils/PartiallySupportedFlowKeywords.cpp

@@ -381,6 +381,12 @@ partiallySupported()
                {13,{false, allow_values<int> {10}, "MESSAGES(PRTGRPMS): option is not supported – will continue"}}, // GROUP_PRINT_LIMIT
             },
          },
+         {
+            "NETBALAN",
+            {
+               {5,{false, allow_values<int> {10}, "NETBALAN(THPMXITE): option is not supported – will continue"}}, // MAX_ITER_THP
+            },
+         },
          {
             "NETWORK",
             {
@@ -579,6 +585,16 @@ partiallySupported()
                {3,{false, allow_values<double> {}, "MULTFLT(FLT-DIF): the diffusivity multiplier option is not supported – will continue"}}, // NOT_DEFINED
             },
          },
+         {
+            "NETBALAN",
+            {
+               {1,{true, [](const double value) { return value <= 0.0; }, "NETBALAN(NSTEP): only negative values or 0 supported – will stop"}}, // TIME_INTERVAL
+               {4,{false, allow_values<double> {0.01}, "NETBALAN(GRPCNV): not supported – will continue"}}, // THP_CONVERGENCE_LIMIT
+               {6,{false, allow_values<double> {1e20}, "NETBALAN(NTRGERR): not supported – will continue"}}, // TARGET_BALANCE_ERROR
+               {7,{false, allow_values<double> {1e20}, "NETBALAN(NMAXERR): not supported – will continue"}}, // MAX_BALANCE_ERROR
+               {8,{false, allow_values<double> {}, "NETBALAN(NTSMIN): not supported – will continue"}}, // MIN_TIME_STEP
+            },
+         },
          {
             "NNC",
             {

opm/simulators/utils/UnsupportedFlowKeywords.cpp

@@ -398,7 +398,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"NARROW", {false, std::nullopt}},
         {"NCONSUMP", {false, std::nullopt}},
         {"NEFAC", {false, std::nullopt}},
-        {"NETBALAN", {false, std::nullopt}},
         {"NETCOMPA", {false, std::nullopt}},
         {"NEXT", {false, std::nullopt}},
         {"NEXTSTPL", {false, std::nullopt}},

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;
-            void 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,
@@ -369,6 +370,10 @@ namespace Opm {
             // and in the well equations.
             void assemble(const int iterationIdx,
                           const double dt);
+            void assembleImpl(const int iterationIdx,
+                              const double dt,
+                              const std::size_t recursion_level,
+                              DeferredLogger& local_deferredLogger);
 
             // called at the end of a time step
             void timeStepSucceeded(const double& simulationTime, const double dt);

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -2102,14 +2102,14 @@ inferLocalShutWells()
     }
 }
 
-void
+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;
+        return { false, 0.0 };
     }
     node_pressures_ = WellGroupHelpers::computeNetworkPressures(network,
                                                                 this->wellState(),
@@ -2119,6 +2119,8 @@ updateNetworkPressures(const int reportStepIdx)
                                                                 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
@@ -2128,10 +2130,19 @@ updateNetworkPressures(const int reportStepIdx)
             if (it != node_pressures_.end()) {
                 // The well belongs to a group with has a network pressure constraint,
                 // set the dynamic THP constraint of the well accordingly.
-                well->setDynamicThpLimit(it->second);
+                const double new_limit = it->second;
+                well->setDynamicThpLimit(new_limit);
+                const 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;
+                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, network_imbalance };
 }
 
 void

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -250,7 +250,7 @@ protected:
 
     bool wasDynamicallyShutThisTimeStep(const int well_index) const;
 
-    void 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

@@ -870,7 +870,6 @@ namespace Opm {
             return;
         }
 
-
         updatePerforationIntensiveQuantities();
 
         if (iterationIdx == 0) {
@@ -886,7 +885,27 @@ namespace Opm {
             OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "assemble() failed (It=0): ",
                                                terminal_output_, grid().comm());
         }
-        updateWellControls(local_deferredLogger, /* check group controls */ true);
+
+        assembleImpl(iterationIdx, dt, 0, local_deferredLogger);
+
+        last_report_.converged = true;
+        last_report_.assemble_time_well += perfTimer.stop();
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    assembleImpl(const int iterationIdx,
+                 const double dt,
+                 const std::size_t recursion_level,
+                 DeferredLogger& local_deferredLogger)
+    {
+
+        const auto [network_changed, network_imbalance] = updateWellControls(local_deferredLogger, /* check group controls */ true);
 
         bool alq_updated = false;
         OPM_BEGIN_PARALLEL_TRY_CATCH();
@@ -911,10 +930,23 @@ namespace Opm {
             WellGroupHelpers::updateGuideRates(fieldGroup, schedule(), summaryState, this->phase_usage_, reportStepIdx, simulationTime,
                                         this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger);
         }
-        last_report_.converged = true;
-        last_report_.assemble_time_well += perfTimer.stop();
+
+        // Maybe do a recursive call to iterate network and well controls.
+        if (network_changed) {
+            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);
+                }
+            }
+        }
     }
 
+
+
+
+
     template<typename TypeTag>
     bool
     BlackoilWellModel<TypeTag>::
@@ -1433,21 +1465,49 @@ namespace Opm {
 
 
     template<typename TypeTag>
-    void
+    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 ;
+        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);
 
-        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;
 
@@ -1500,6 +1560,8 @@ namespace Opm {
         // update wsolvent fraction for REIN wells
         const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
         updateWsolvent(fieldGroup, episodeIdx,  this->nupcolWellState());
+
+        return { network_changed, network_imbalance };
     }
 
 

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -347,6 +347,11 @@ void WellInterfaceGeneric::setDynamicThpLimit(const double thp_limit)
     dynamic_thp_limit_ = thp_limit;
 }
 
+std::optional<double> WellInterfaceGeneric::getDynamicThpLimit() const
+{
+    return dynamic_thp_limit_;
+}
+
 void WellInterfaceGeneric::updatePerforatedCell(std::vector<bool>& is_cell_perforated)
 {
 

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -101,6 +101,7 @@ public:
     void setRepRadiusPerfLength();
     void setWsolvent(const double wsolvent);
     void setDynamicThpLimit(const double thp_limit);
+    std::optional<double> getDynamicThpLimit() const;
     void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
 
     /// Returns true if the well has one or more THP limits/constraints.