Add network model iterations.

Implemented by conditionally recursing in BlackoilWellModel::assembleImpl().

opm/simulators/wells/BlackoilWellModel.hpp

@@ -282,7 +282,7 @@ 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 updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls);
 
             void updateAndCommunicate(const int reportStepIdx,
                                       const int iterationIdx,
@@ -369,6 +369,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 int 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

@@ -2092,14 +2092,14 @@ inferLocalShutWells()
     }
 }
 
-void
+bool
 BlackoilWellModelGeneric::
 updateNetworkPressures(const int reportStepIdx)
 {
     // Get the network and return if inactive.
     const auto& network = schedule()[reportStepIdx].network();
     if (!network.active()) {
-        return;
+        return false;
     }
     node_pressures_ = WellGroupHelpers::computeNetworkPressures(network,
                                                                 this->wellState(),
@@ -2109,6 +2109,7 @@ updateNetworkPressures(const int reportStepIdx)
                                                                 reportStepIdx);
 
     // Set the thp limits of wells
+    bool active_limit_change = false;
     for (auto& well : well_container_generic_) {
         // Producers only, since we so far only support the
         // "extended" network model (properties defined by
@@ -2118,10 +2119,18 @@ 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;
+                }
             }
         }
     }
+    return active_limit_change;
 }
 
 void

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -250,7 +250,7 @@ protected:
 
     bool wasDynamicallyShutThisTimeStep(const int well_index) const;
 
-    void updateNetworkPressures(const int reportStepIdx);
+    bool 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 int recursion_level,
+                 DeferredLogger& local_deferredLogger)
+    {
+
+        const bool network_changed = 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) {
+            const int max_local_network_iterations = 3;
+            const int last_newton_iteration_to_iterate_network = 2;
+            if (iterationIdx <= last_newton_iteration_to_iterate_network) {
+                if (recursion_level < max_local_network_iterations) {
+                    assembleImpl(iterationIdx, dt, recursion_level + 1, local_deferredLogger);
+                }
+            }
+        }
     }
 
+
+
+
+
     template<typename TypeTag>
     bool
     BlackoilWellModel<TypeTag>::
@@ -1433,21 +1465,21 @@ namespace Opm {
 
 
     template<typename TypeTag>
-    void
+    bool
     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;
 
         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 bool network_changed = updateNetworkPressures(episodeIdx);
 
         std::set<std::string> switched_wells;
 
@@ -1500,6 +1532,8 @@ namespace Opm {
         // update wsolvent fraction for REIN wells
         const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
         updateWsolvent(fieldGroup, episodeIdx,  this->nupcolWellState());
+
+        return network_changed;
     }