Merge pull request #4598 from GitPaean/network_continue

various network improvements

opm/simulators/wells/BlackoilWellModel.hpp

@@ -293,7 +293,8 @@ namespace Opm {
             void prepareTimeStep(DeferredLogger& deferred_logger);
             void initPrimaryVariablesEvaluation() const;
 
-            std::pair<bool, bool> updateWellControls(DeferredLogger& deferred_logger, const std::size_t network_update_it = 0);
+            std::pair<bool, bool>
+            updateWellControls(const bool mandatory_network_balance, DeferredLogger& deferred_logger, const bool relax_network_tolerance = false);
 
             void updateAndCommunicate(const int reportStepIdx,
                                       const int iterationIdx,
@@ -367,6 +368,10 @@ namespace Opm {
 
             SimulatorReportSingle last_report_{};
 
+            // solve to get a good network solution, group and well states might be updated during the process.
+            // the reservoir should stay static during this solution procedure.
+            void balanceNetwork(DeferredLogger& deferred_logger);
+
             // used to better efficiency of calcuation
             mutable BVector scaleAddRes_{};
 
@@ -390,11 +395,14 @@ namespace Opm {
             // 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,
+            std::pair<bool, bool> updateWellControlsAndNetworkIteration(const bool mandatory_network_balance,
+                                                                        const bool relax_network_tolerance,
+                                                                        const double dt,
                                                                         DeferredLogger& local_deferredLogger);
 
-            bool updateWellControlsAndNetwork(const double dt, DeferredLogger& local_deferredLogger);
+            bool updateWellControlsAndNetwork(const bool mandatory_network_balance,
+                                              const double dt,
+                                              DeferredLogger& local_deferredLogger);
 
             // called at the end of a time step
             void timeStepSucceeded(const double& simulationTime, const double dt);

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>
@@ -905,6 +906,30 @@ hasTHPConstraints() const
     return BlackoilWellModelConstraints(*this).hasTHPConstraints();
 }
 
+bool
+BlackoilWellModelGeneric::
+needRebalanceNetwork(const int report_step) const
+{
+    const auto& network = schedule()[report_step].network();
+    if (!network.active()) {
+        return false;
+    }
+
+    bool network_rebalance_necessary = false;
+    for (const auto& well : well_container_generic_) {
+        const auto& events = this->wellState().well(well->indexOfWell()).events;
+        const bool is_partof_network = network.has_node(well->wellEcl().groupName());
+        // TODO: we might find more relevant events to be included here
+        if (is_partof_network && events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE)) {
+            network_rebalance_necessary = true;
+            break;
+        }
+    }
+    network_rebalance_necessary = comm_.max(network_rebalance_necessary);
+
+    return network_rebalance_necessary;
+}
+
 bool
 BlackoilWellModelGeneric::
 forceShutWellByName(const std::string& wellname,
@@ -958,15 +983,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 +1002,40 @@ updateNetworkPressures(const int reportStepIdx)
                                                                 schedule(),
                                                                 reportStepIdx);
 
-    // Set the thp limits of wells
-    bool active_limit_change = false;
-    double network_imbalance = 0.0;
+    // here, the network imbalance is the difference between the previous nodal pressure and the new nodal pressure
+    double network_imbalance = 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.
@@ -988,17 +1046,16 @@ updateNetworkPressures(const int reportStepIdx)
                 // set the dynamic THP constraint of the well accordingly.
                 const double new_limit = it->second;
                 well->setDynamicThpLimit(new_limit);
-                const SingleWellState& ws = this->wellState()[well->indexOfWell()];
+                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));
+                // TODO: not sure why the thp is NOT updated properly elsewhere
+                if (thp_is_limit) {
+                    ws.thp = well->getTHPConstraint(summaryState_);
                 }
             }
         }
     }
-    return { active_limit_change, network_imbalance };
+    return network_imbalance;
 }
 
 void
@@ -1238,6 +1295,12 @@ bool
 BlackoilWellModelGeneric::
 shouldBalanceNetwork(const int reportStepIdx, const int iterationIdx) const
 {
+    // if network is not active, we do not need to balance the network
+    const auto& network = schedule()[reportStepIdx].network();
+    if (!network.active()) {
+        return false;
+    }
+
     const auto& balance = schedule()[reportStepIdx].network_balance();
     if (balance.mode() == Network::Balance::CalcMode::TimeStepStart) {
         return iterationIdx == 0;
@@ -1253,17 +1316,7 @@ shouldBalanceNetwork(const int reportStepIdx, const int iterationIdx) const
     }
 }
 
-bool
-BlackoilWellModelGeneric::
-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 iteration < balance.pressure_max_iter() &&
-           network_imbalance > balance.pressure_tolerance();
-}
+
 
 std::vector<int>
 BlackoilWellModelGeneric::

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -158,6 +158,9 @@ public:
     /// Return true if any well has a THP constraint.
     bool hasTHPConstraints() const;
 
+    /// Whether it is necessary to re-balance network
+    bool needRebalanceNetwork(const int report_step) const;
+
     /// Shut down any single well
     /// Returns true if the well was actually found and shut.
     bool forceShutWellByName(const std::string& wellname,
@@ -177,10 +180,6 @@ public:
     bool shouldBalanceNetwork(const int reportStepIndex,
                               const int iterationIdx) const;
 
-    bool moreNetworkIteration(const int reportStepIdx,
-                              const std::size_t iteration,
-                              const double network_imbalance) const;
-
     template<class Serializer>
     void serializeOp(Serializer& serializer)
     {
@@ -302,7 +301,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

@@ -25,12 +25,13 @@
 #include <opm/grid/utility/cartesianToCompressed.hpp>
 
 #include <opm/input/eclipse/Units/UnitSystem.hpp>
+#include <opm/input/eclipse/Schedule/Network/Balance.hpp>
+#include <opm/input/eclipse/Schedule/Network/ExtNetwork.hpp>
 
 #include <opm/simulators/wells/BlackoilWellModelConstraints.hpp>
 #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>
@@ -746,7 +747,25 @@ namespace Opm {
 
         well_container_generic_.clear();
         for (auto& w : well_container_)
-          well_container_generic_.push_back(w.get());
+            well_container_generic_.push_back(w.get());
+
+        const auto& network = schedule()[time_step].network();
+        if (network.active() && !this->node_pressures_.empty()) {
+            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.
+                if (well->isProducer()) {
+                    const auto it = node_pressures_.find(well->wellEcl().groupName());
+                    if (it != node_pressures_.end()) {
+                        // The well belongs to a group which has a network nodal pressure,
+                        // set the dynamic THP constraint based on the network nodal pressure
+                        const double nodal_pressure = it->second;
+                        well->setDynamicThpLimit(nodal_pressure);
+                    }
+                }
+            }
+        }
     }
 
 
@@ -831,6 +850,43 @@ namespace Opm {
 
 
 
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    balanceNetwork(DeferredLogger& deferred_logger) {
+        const double dt = this->ebosSimulator_.timeStepSize();
+        // TODO: should we also have the group and network backed-up here in case the solution did not get converged?
+        auto& well_state = this->wellState();
+        constexpr size_t max_iter = 100;
+        bool converged = false;
+        size_t iter = 0;
+        bool changed_well_group = false;
+        do {
+            changed_well_group = updateWellControlsAndNetwork(true, dt, deferred_logger);
+            assembleWellEqWithoutIteration(dt, deferred_logger);
+            converged = this->getWellConvergence(this->B_avg_, true).converged() && !changed_well_group;
+            if (converged) {
+                break;
+            }
+            ++iter;
+            for (auto& well : this->well_container_) {
+                const auto& summary_state = this->ebosSimulator_.vanguard().summaryState();
+                well->solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
+            }
+            this->initPrimaryVariablesEvaluation();
+        } while (iter < max_iter);
+
+        if (!converged) {
+            const std::string msg = fmt::format("balanceNetwork did not get converged with {} iterations, and unconverged "
+                                                "network balance result will be used", max_iter);
+            deferred_logger.warning(msg);
+        } else {
+            const std::string msg = fmt::format("balanceNetwork get converged with {} iterations", iter);
+            deferred_logger.debug(msg);
+        }
+    }
+
+
 
 
     template<typename TypeTag>
@@ -850,13 +906,17 @@ namespace Opm {
         Dune::Timer perfTimer;
         perfTimer.start();
 
-        if ( ! wellsActive() ) {
-            return;
+        {
+            const int episodeIdx = ebosSimulator_.episodeIndex();
+            const auto& network = schedule()[episodeIdx].network();
+            if ( !wellsActive() && !network.active() ) {
+                return;
+            }
         }
 
         updatePerforationIntensiveQuantities();
 
-        if (iterationIdx == 0) {
+        if (iterationIdx == 0 && wellsActive()) {
             // try-catch is needed here as updateWellControls
             // contains global communication and has either to
             // be reached by all processes or all need to abort
@@ -871,7 +931,13 @@ namespace Opm {
                                            terminal_output_, grid().comm());
         }
 
-        const bool well_group_control_changed = updateWellControlsAndNetwork(dt, local_deferredLogger);
+        const bool well_group_control_changed = updateWellControlsAndNetwork(false, dt, local_deferredLogger);
+
+        // even when there is no wells active, the network nodal pressure still need to be updated through updateWellControlsAndNetwork()
+        // but there is no need to assemble the well equations
+        if ( ! wellsActive() ) {
+            return;
+        }
 
         assembleWellEqWithoutIteration(dt, local_deferredLogger);
 
@@ -887,16 +953,29 @@ namespace Opm {
     template<typename TypeTag>
     bool
     BlackoilWellModel<TypeTag>::
-    updateWellControlsAndNetwork(const double dt, DeferredLogger& local_deferredLogger)
+    updateWellControlsAndNetwork(const bool mandatory_network_balance, 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;
+        // after certain number of the iterations, we use relaxed tolerance for the network update
+        constexpr size_t iteration_to_relax = 100;
+        // after certain number of the iterations, we terminate
+        constexpr size_t max_iteration = 200;
         std::size_t network_update_iteration = 0;
         while (do_network_update) {
+            if (network_update_iteration == iteration_to_relax) {
+                local_deferredLogger.info(" we begin using relaxed tolerance for network update now after " + std::to_string(iteration_to_relax) + " iterations ");
+            }
+            const bool relax_network_balance = network_update_iteration >= iteration_to_relax;
             std::tie(do_network_update, well_group_control_changed) =
-                    updateWellControlsAndNetworkIteration(dt, network_update_iteration, local_deferredLogger);
+                    updateWellControlsAndNetworkIteration(mandatory_network_balance, relax_network_balance, dt,local_deferredLogger);
             ++network_update_iteration;
+            if (network_update_iteration >= max_iteration) {
+                local_deferredLogger.info("maximum of " + std::to_string(max_iteration) + " iterations has been used, we stop the network update now, "
+                                                                                          "the simulation will continue with unconvergeed network results");
+                break;
+            }
         }
         return well_group_control_changed;
     }
@@ -907,11 +986,13 @@ namespace Opm {
     template<typename TypeTag>
     std::pair<bool, bool>
     BlackoilWellModel<TypeTag>::
-    updateWellControlsAndNetworkIteration(const double dt,
-                                          const std::size_t network_update_iteration,
+    updateWellControlsAndNetworkIteration(const bool mandatory_network_balance,
+                                          const bool relax_network_tolerance,
+                                          const double dt,
                                           DeferredLogger& local_deferredLogger)
     {
-        auto [well_group_control_changed, more_network_update] = updateWellControls(local_deferredLogger, network_update_iteration);
+        auto [well_group_control_changed, more_network_update] =
+                updateWellControls(mandatory_network_balance, local_deferredLogger, relax_network_tolerance);
 
         bool alq_updated = false;
         OPM_BEGIN_PARALLEL_TRY_CATCH();
@@ -1490,28 +1571,27 @@ namespace Opm {
     template<typename TypeTag>
     std::pair<bool, bool>
     BlackoilWellModel<TypeTag>::
-    updateWellControls(DeferredLogger& deferred_logger,
-                       const std::size_t network_update_it)
+    updateWellControls(const bool mandatory_network_balance, DeferredLogger& deferred_logger, const bool relax_network_tolerance)
     {
-        // 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 };
-
         const int episodeIdx = ebosSimulator_.episodeIndex();
+        const auto& network = schedule()[episodeIdx].network();
+        if (!wellsActive() && !network.active()) {
+            return {false, false};
+        }
+
         const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
         const auto& comm = ebosSimulator_.vanguard().grid().comm();
         updateAndCommunicateGroupData(episodeIdx, iterationIdx);
 
         // 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);
+        if (shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) {
+            const auto local_network_imbalance = updateNetworkPressures(episodeIdx);
             const double network_imbalance = comm.max(local_network_imbalance);
-            if (network_changed) {
-                more_network_update = moreNetworkIteration(episodeIdx, network_update_it, network_imbalance);
-            }
+            const auto& balance = schedule()[episodeIdx].network_balance();
+            constexpr double relaxtion_factor = 10.0;
+            const double tolerance = relax_network_tolerance ? relaxtion_factor * balance.pressure_tolerance() : balance.pressure_tolerance();
+            more_network_update = network_imbalance > tolerance;
         }
 
         bool changed_well_group = false;
@@ -1724,13 +1804,13 @@ namespace Opm {
 
 
 
-
-
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
     prepareTimeStep(DeferredLogger& deferred_logger)
     {
+        const bool network_rebalance_necessary = this->needRebalanceNetwork(ebosSimulator_.episodeIndex());
+
         for (const auto& well : well_container_) {
             auto& events = this->wellState().well(well->indexOfWell()).events;
             if (events.hasEvent(WellState::event_mask)) {
@@ -1753,6 +1833,11 @@ namespace Opm {
             }
         }
         updatePrimaryVariables(deferred_logger);
+
+        if (network_rebalance_necessary) {
+            // this is to obtain good network solution
+            balanceNetwork(deferred_logger);
+        }
     }
 
     template<typename TypeTag>

opm/simulators/wells/MultisegmentWell.hpp

@@ -124,6 +124,10 @@ namespace Opm
                                     const WellState& well_state,
                                     DeferredLogger& deferred_logger) override;
 
+        virtual void solveEqAndUpdateWellState(const SummaryState& summary_state,
+                                               WellState& well_state,
+                                               DeferredLogger& deferred_logger) override; // const?
+
         virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
                                                  const WellState& well_state,
                                                  DeferredLogger& deferred_logger) override; // should be const?

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -531,6 +531,26 @@ namespace Opm
 
 
 
+    template <typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    solveEqAndUpdateWellState(const SummaryState& summary_state,
+                              WellState& well_state,
+                              DeferredLogger& deferred_logger)
+    {
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
+
+        // We assemble the well equations, then we check the convergence,
+        // which is why we do not put the assembleWellEq here.
+        const BVectorWell dx_well = this->linSys_.solve();
+
+        updateWellState(summary_state, dx_well, well_state, deferred_logger);
+    }
+
+
+
+
+
     template <typename TypeTag>
     void
     MultisegmentWell<TypeTag>::

opm/simulators/wells/StandardWell.hpp

@@ -172,9 +172,9 @@ namespace Opm
                                     const WellState& well_state,
                                     DeferredLogger& deferred_logger) override;
 
-        void solveEqAndUpdateWellState(const SummaryState& summary_state,
-                                       WellState& well_state,
-                                       DeferredLogger& deferred_logger);
+        virtual void solveEqAndUpdateWellState(const SummaryState& summary_state,
+                                               WellState& well_state,
+                                               DeferredLogger& deferred_logger) override;
 
         virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
                                                  const WellState& well_state,

opm/simulators/wells/WellInterface.hpp

@@ -160,6 +160,10 @@ public:
                                                  DeferredLogger& deferred_logger,
                                                  const bool relax_tolerance) const = 0;
 
+    virtual void solveEqAndUpdateWellState(const SummaryState& summary_state,
+                                           WellState& well_state,
+                                           DeferredLogger& deferred_logger) = 0;
+
     void assembleWellEq(const Simulator& ebosSimulator,
                         const double dt,
                         WellState& well_state,

parallelTests.cmake

@@ -141,6 +141,14 @@ add_test_compare_parallel_simulation(CASENAME aquflux_01
                                      DIR aquifers
                                      TEST_ARGS --enable-tuning=true)
 
+add_test_compare_parallel_simulation(CASENAME network_balance_01
+		                             FILENAME NETWORK-01
+		                             SIMULATOR flow
+		                             ABS_TOL 0.04
+		                             REL_TOL 0.02
+		                             DIR network
+		                             TEST_ARGS --enable-tuning=true)
+
 add_test_compare_parallel_simulation(CASENAME numerical_aquifer_3d_1aqu
                                      FILENAME 3D_1AQU_3CELLS
                                      SIMULATOR flow

regressionTests.cmake

@@ -92,6 +92,14 @@ add_test_compareECLFiles(CASENAME spe1_precsalt
                          REL_TOL ${rel_tol}
                          DIR spe1_precsalt)
 
+add_test_compareECLFiles(CASENAME network_balance_01
+		                 FILENAME NETWORK-01
+		                 SIMULATOR flow
+		                 ABS_TOL ${abs_tol}
+		                 REL_TOL ${rel_tol}
+		                 DIR network
+		                 TEST_ARGS --enable-tuning=true)
+
 add_test_compareECLFiles(CASENAME gas_precsalt
                          FILENAME GASWATER_VAPWAT_PRECSALT
                          SIMULATOR flow

tests/test_glift1.cpp

@@ -147,7 +147,7 @@ BOOST_AUTO_TEST_CASE(G1)
     Opm::DeferredLogger deferred_logger;
     well_model.calculateExplicitQuantities(deferred_logger);
     well_model.prepareTimeStep(deferred_logger);
-    well_model.updateWellControls(deferred_logger);
+    well_model.updateWellControls(false, deferred_logger);
     well_model.initPrimaryVariablesEvaluation();
     Opm::WellInterface<TypeTag> *well_ptr = well_model.getWell("B-1H").get();
     StdWell *std_well = dynamic_cast<StdWell *>(well_ptr);

tests/update_reference_data.sh

@@ -77,6 +77,7 @@ tests[spe1_water]="flow spe1 SPE1CASE1_WATER"
 tests[spe1_spider]="flow radial_grid SPIDER_CAKESLICE"
 tests[spe1_radial]="flow radial_grid RADIAL_CAKESLICE"
 tests[spe1_import]="flow spe1 SPE1CASE1_IMPORT"
+tests[network_balance_01]="flow network NETWORK-01"
 tests[jfunc_01]="flow jfunc JFUNC-01"
 tests[pinch_nopinch_1x1x10]="flow pinch PINCH10_NOPINCH"
 tests[ctaquifer_2d_oilwater]="flow aquifer-oilwater 2D_OW_CTAQUIFER"