Merge pull request #4596 from GitPaean/refactor_network_update_test

refactor the network update

opm/simulators/wells/BlackoilWellModel.hpp

@@ -292,7 +292,8 @@ namespace Opm {
             // at the beginning of each time step (Not report step)
             void prepareTimeStep(DeferredLogger& deferred_logger);
             void initPrimaryVariablesEvaluation() const;
-            std::tuple<bool, bool, double> updateWellControls(DeferredLogger& deferred_logger);
+
+            std::pair<bool, bool> updateWellControls(DeferredLogger& deferred_logger, const std::size_t network_update_it = 0);
 
             void updateAndCommunicate(const int reportStepIdx,
                                       const int iterationIdx,
@@ -384,10 +385,16 @@ namespace Opm {
             // and in the well equations.
             void assemble(const int iterationIdx,
                           const double dt);
-            bool assembleImpl(const int iterationIdx,
-                              const double dt,
-                              const std::size_t recursion_level,
-                              DeferredLogger& local_deferredLogger);
+
+            // well controls and network pressures affect each other and are solved in an iterative manner.
+            // 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 double dt,
+                                                                        const std::size_t network_update_iteration,
+                                                                        DeferredLogger& local_deferredLogger);
+
+            bool updateWellControlsAndNetwork(const double dt, DeferredLogger& local_deferredLogger);
 
             // called at the end of a time step
             void timeStepSucceeded(const double& simulationTime, const double dt);
@@ -424,6 +431,11 @@ namespace Opm {
 
             void assembleWellEq(const double dt, DeferredLogger& deferred_logger);
 
+            void prepareWellsBeforeAssembling(const double dt, DeferredLogger& deferred_logger);
+
+            // TODO: finding a better naming
+            void assembleWellEqWithoutIteration(const double dt, DeferredLogger& deferred_logger);
+
             bool maybeDoGasLiftOptimize(DeferredLogger& deferred_logger);
 
             void gasLiftOptimizationStage1(DeferredLogger& deferred_logger,

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1255,13 +1255,13 @@ shouldBalanceNetwork(const int reportStepIdx, const int iterationIdx) const
 
 bool
 BlackoilWellModelGeneric::
-shouldIterateNetwork(const int reportStepIdx,
-                     const std::size_t recursion_level,
+moreNetworkIteration(const int reportStepIdx,
+                     const std::size_t iteration,
                      const double network_imbalance) const
 {
     const auto& balance = schedule()[reportStepIdx].network_balance();
     // Iterate if not converged, and number of iterations is not yet max (NETBALAN item 3).
-    return recursion_level < balance.pressure_max_iter() &&
+    return iteration < balance.pressure_max_iter() &&
            network_imbalance > balance.pressure_tolerance();
 }
 

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -177,8 +177,8 @@ public:
     bool shouldBalanceNetwork(const int reportStepIndex,
                               const int iterationIdx) const;
 
-    bool shouldIterateNetwork(const int reportStepIndex,
-                              const std::size_t recursion_level,
+    bool moreNetworkIteration(const int reportStepIdx,
+                              const std::size_t iteration,
                               const double network_imbalance) const;
 
     template<class Serializer>

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -30,6 +30,7 @@
 #include <opm/simulators/wells/VFPProperties.hpp>
 #include <opm/simulators/utils/MPIPacker.hpp>
 #include <opm/simulators/linalg/bda/WellContributions.hpp>
+#include <opm/input/eclipse/Schedule/Network/Balance.hpp>
 
 #if HAVE_MPI
 #include <ebos/eclmpiserializer.hh>
@@ -870,7 +871,9 @@ namespace Opm {
                                            terminal_output_, grid().comm());
         }
 
-        const bool well_group_control_changed = assembleImpl(iterationIdx, dt, 0, local_deferredLogger);
+        const bool well_group_control_changed = updateWellControlsAndNetwork(dt, local_deferredLogger);
+
+        assembleWellEqWithoutIteration(dt, local_deferredLogger);
 
         // if group or well control changes we don't consider the
         // case converged
@@ -881,17 +884,34 @@ namespace Opm {
 
 
 
-
     template<typename TypeTag>
     bool
     BlackoilWellModel<TypeTag>::
-    assembleImpl(const int iterationIdx,
-                 const double dt,
-                 const std::size_t recursion_level,
-                 DeferredLogger& local_deferredLogger)
+    updateWellControlsAndNetwork(const double dt, DeferredLogger& local_deferredLogger)
     {
+        // not necessarily that we always need to update once of the network solutions
+        bool do_network_update = true;
+        bool well_group_control_changed = false;
+        std::size_t network_update_iteration = 0;
+        while (do_network_update) {
+            std::tie(do_network_update, well_group_control_changed) =
+                    updateWellControlsAndNetworkIteration(dt, network_update_iteration, local_deferredLogger);
+            ++network_update_iteration;
+        }
+        return well_group_control_changed;
+    }
 
-        auto [well_group_control_changed, network_changed, network_imbalance] = updateWellControls(local_deferredLogger);
+
+
+
+    template<typename TypeTag>
+    std::pair<bool, bool>
+    BlackoilWellModel<TypeTag>::
+    updateWellControlsAndNetworkIteration(const double dt,
+                                          const std::size_t network_update_iteration,
+                                          DeferredLogger& local_deferredLogger)
+    {
+        auto [well_group_control_changed, more_network_update] = updateWellControls(local_deferredLogger, network_update_iteration);
 
         bool alq_updated = false;
         OPM_BEGIN_PARALLEL_TRY_CATCH();
@@ -900,9 +920,10 @@ namespace Opm {
             initPrimaryVariablesEvaluation();
 
             alq_updated = maybeDoGasLiftOptimize(local_deferredLogger);
-            assembleWellEq(dt, local_deferredLogger);
+
+            prepareWellsBeforeAssembling(dt, local_deferredLogger);
         }
-        OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "assemble() failed: ",
+        OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "updateWellControlsAndNetworkIteration() failed: ",
                                        terminal_output_, grid().comm());
 
         //update guide rates
@@ -915,23 +936,16 @@ namespace Opm {
             std::vector<double> pot(numPhases(), 0.0);
             const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
             WellGroupHelpers::updateGuideRates(fieldGroup, schedule(), summaryState, this->phase_usage_, reportStepIdx, simulationTime,
-                                        this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger);
+                                               this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger);
         }
 
-        // Maybe do a recursive call to iterate network and well controls.
-        if (network_changed) {
-            if (shouldBalanceNetwork(reportStepIdx, iterationIdx) &&
-                shouldIterateNetwork(reportStepIdx, recursion_level, network_imbalance)) {
-                well_group_control_changed = assembleImpl(iterationIdx, dt, recursion_level + 1, local_deferredLogger);
-            }
-        }
-        return well_group_control_changed;
+
+        return {more_network_update, well_group_control_changed};
     }
 
 
 
 
-
     template<typename TypeTag>
     bool
     BlackoilWellModel<TypeTag>::
@@ -1143,6 +1157,29 @@ namespace Opm {
         }
     }
 
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    prepareWellsBeforeAssembling(const double dt, DeferredLogger& deferred_logger)
+    {
+        for (auto& well : well_container_) {
+            well->prepareWellBeforeAssembling(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger);
+        }
+    }
+
+
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    assembleWellEqWithoutIteration(const double dt, DeferredLogger& deferred_logger)
+    {
+        for (auto& well: well_container_) {
+            well->assembleWellEqWithoutIteration(ebosSimulator_, dt, this->wellState(), this->groupState(),
+                                                 deferred_logger);
+        }
+    }
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
@@ -1451,25 +1488,31 @@ namespace Opm {
 
 
     template<typename TypeTag>
-    std::tuple<bool, bool, double>
+    std::pair<bool, bool>
     BlackoilWellModel<TypeTag>::
-    updateWellControls(DeferredLogger& deferred_logger)
+    updateWellControls(DeferredLogger& deferred_logger,
+                       const std::size_t network_update_it)
     {
         // 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, false, 0.0 };
+        if( !wellsActive() ) return { false, false };
 
         const int episodeIdx = ebosSimulator_.episodeIndex();
         const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
         const auto& comm = ebosSimulator_.vanguard().grid().comm();
         updateAndCommunicateGroupData(episodeIdx, iterationIdx);
 
-        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);
+        // network related
+        bool more_network_update = false;
+        if (shouldBalanceNetwork(episodeIdx, iterationIdx)) {
+            const auto [local_network_changed, 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);
+            }
+        }
 
         bool changed_well_group = false;
         // Check group individual constraints.
@@ -1515,7 +1558,7 @@ namespace Opm {
         const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
         updateWsolvent(fieldGroup, episodeIdx,  this->nupcolWellState());
 
-        return { changed_well_group, network_changed, network_imbalance };
+        return { changed_well_group, more_network_update };
     }
 
 

opm/simulators/wells/WellInterface.hpp

@@ -166,6 +166,20 @@ public:
                         const GroupState& group_state,
                         DeferredLogger& deferred_logger);
 
+    void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
+                                        const double dt,
+                                        WellState& well_state,
+                                        const GroupState& group_state,
+                                        DeferredLogger& deferred_logger);
+
+    // TODO: better name or further refactoring the function to make it more clear
+    void prepareWellBeforeAssembling(const Simulator& ebosSimulator,
+                                     const double dt,
+                                     WellState& well_state,
+                                     const GroupState& group_state,
+                                     DeferredLogger& deferred_logger);
+
+
     virtual void computeWellRatesWithBhp(
         const Simulator& ebosSimulator,
         const double& bhp,

opm/simulators/wells/WellInterface_impl.hpp

@@ -470,6 +470,42 @@ namespace Opm
                    const GroupState& group_state,
                    DeferredLogger& deferred_logger)
     {
+
+        prepareWellBeforeAssembling(ebosSimulator, dt, well_state, group_state, deferred_logger);
+
+        assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, group_state, deferred_logger);
+    }
+
+
+
+    template <typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
+                                   const double dt,
+                                   WellState& well_state,
+                                   const GroupState& group_state,
+                                   DeferredLogger& deferred_logger)
+    {
+        const auto& summary_state = ebosSimulator.vanguard().summaryState();
+        const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
+        const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
+        // TODO: the reason to have inj_controls and prod_controls in the arguments, is that we want to change the control used for the well functions
+        // TODO: maybe we can use std::optional or pointers to simplify here
+        assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
+    }
+
+
+
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    prepareWellBeforeAssembling(const Simulator& ebosSimulator,
+                                const double dt,
+                                WellState& well_state,
+                                const GroupState& group_state,
+                                DeferredLogger& deferred_logger)
+    {
         const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
 
         if (param_.check_well_operability_iter_)
@@ -522,11 +558,6 @@ namespace Opm
             changed_to_stopped_this_step_ = false;
             this->changed_to_open_this_step_ = true;
         }
-
-        const auto& summary_state = ebosSimulator.vanguard().summaryState();
-        const auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
-        const auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
-        assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
     }
 
     template<typename TypeTag>